Alpha diversity
label_map <- c(
"Control" = "Cold-control",
"Treatment" = "Cold-intervention",
"richness" = "Species Richness",
"neutral" = "Neutral Diversity",
"phylogenetic" = "Phylogenetic Diversity"
)
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point %in% c("FMT1","Acclimation", "FMT2")) %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5, show.legend = FALSE) +
facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
# facet_nested(.~metric+type, labeller = labeller(metric = label_map, type = label_map))+
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
theme_minimal() +
theme(axis.text.x=element_text(size=10))+
labs(x = "Time Point", y = "value")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT1"|time_point=="Acclimation"|time_point=="FMT2")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(5.7695) ( log )
Formula: richness ~ type * time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
757.7 783.7 -367.8 735.7 68
=======
load("data/data_27022025.Rdata")
load("data/calculations_28022025.Rdata")
What is the effect of FMT and temperature on the GM after 1 week compared to acclimation?
CI vs CC
Alpha diversity
label_map <- c(
"Control" = "Cold-control",
"Treatment" = "Cold-intervention",
"richness" = "Species Richness",
"neutral" = "Neutral Diversity",
"phylogenetic" = "Phylogenetic Diversity"
)
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(type %in% c("Control","Treatment") & time_point %in% c("FMT1","Acclimation")) %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5, show.legend = FALSE) +
facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
# facet_nested(.~metric+type, labeller = labeller(metric = label_map, type = label_map))+
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
values=c("#4477AA", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control", "Treatment"),
values=c("#4477AA50","#76b18350")) +
theme_minimal() +
theme(axis.text.x=element_text(size=10))+
labs(x = "Time Point", y = "value")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Hot_control")
Richness
Modelq0GLMMNB <- MASS::glm.nb(richness ~ type*time_point, data = alpha_div_meta,trace=TRUE)
summary(Modelq0GLMMNB)
Call:
MASS::glm.nb(formula = richness ~ type * time_point, data = alpha_div_meta,
trace = TRUE, init.theta = 2.532718573, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.9512 0.2275 17.368 <2e-16 ***
typeTreatment -0.1372 0.3132 -0.438 0.661
time_pointFMT1 -0.3225 0.3140 -1.027 0.304
typeTreatment:time_pointFMT1 0.4500 0.4435 1.015 0.310
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(2.5327) family taken to be 1)
Null deviance: 37.769 on 33 degrees of freedom
Residual deviance: 36.407 on 30 degrees of freedom
AIC: 327.2
Number of Fisher Scoring iterations: 1
Theta: 2.533
Std. Err.: 0.630
2 x log-likelihood: -317.205
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ type)
$emmeans
type emmean SE df asymp.LCL asymp.UCL
Control 3.79 0.157 Inf 3.48 4.10
Treatment 3.88 0.157 Inf 3.57 4.18
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Control - Treatment -0.0878 0.222 Inf -0.396 0.6921
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Neutral
Modelq1n <- lm(formula = neutral ~ type*time_point, data = alpha_div_meta)
summary(Modelq1n)
Call:
lm(formula = neutral ~ type * time_point, data = alpha_div_meta)
Residuals:
Min 1Q Median 3Q Max
-20.699 -6.754 1.842 6.297 15.828
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 21.373 3.716 5.751 2.8e-06 ***
typeTreatment -3.974 5.107 -0.778 0.443
time_pointFMT1 -4.014 5.107 -0.786 0.438
typeTreatment:time_pointFMT1 11.478 7.223 1.589 0.123
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.51 on 30 degrees of freedom
Multiple R-squared: 0.08999, Adjusted R-squared: -0.001011
F-statistic: 0.9889 on 3 and 30 DF, p-value: 0.4113
Phylogenetic
Model_phylo <- lm(formula = phylogenetic ~ type*time_point, data = alpha_div_meta)
summary(Model_phylo)
Call:
lm(formula = phylogenetic ~ type * time_point, data = alpha_div_meta)
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
Scaled residuals:
Min 1Q Median 3Q Max
-2.01194 -0.49654 0.06198 0.55787 1.90952
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0.02586 0.1608
Number of obs: 79, groups: individual, 27
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.9089 0.1696 23.046 <2e-16 ***
typeHot_control 0.5568 0.2283 2.439 0.0147 *
typeTreatment -0.1233 0.2297 -0.537 0.5916
time_pointFMT1 -0.3140 0.2167 -1.449 0.1473
time_pointFMT2 0.1487 0.2173 0.684 0.4937
typeHot_control:time_pointFMT1 -0.1155 0.2983 -0.387 0.6987
typeTreatment:time_pointFMT1 0.4309 0.3068 1.404 0.1602
typeHot_control:time_pointFMT2 -0.3876 0.2988 -1.297 0.1945
typeTreatment:time_pointFMT2 0.4093 0.2990 1.369 0.1710
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) typHt_ typTrt t_FMT1 t_FMT2 tH_:_FMT1 tT:_FMT1 tH_:_FMT2
typHt_cntrl -0.741
typeTretmnt -0.721 0.535
tim_pntFMT1 -0.677 0.503 0.497
tim_pntFMT2 -0.702 0.520 0.507 0.532
typH_:_FMT1 0.496 -0.671 -0.362 -0.727 -0.389
typTr:_FMT1 0.482 -0.357 -0.664 -0.707 -0.378 0.515
typH_:_FMT2 0.515 -0.684 -0.370 -0.388 -0.730 0.516 0.276
typTr:_FMT2 0.498 -0.370 -0.686 -0.384 -0.719 0.280 0.511 0.524
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ type)
$emmeans
type emmean SE df asymp.LCL asymp.UCL
Control 3.85 0.1050 Inf 3.65 4.06
Hot_control 4.24 0.0996 Inf 4.05 4.44
Treatment 4.01 0.1030 Inf 3.81 4.21
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Control - Hot_control -0.389 0.144 Inf -2.711 0.0184
Control - Treatment -0.157 0.145 Inf -1.083 0.5246
Hot_control - Treatment 0.232 0.143 Inf 1.627 0.2342
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 3 estimates
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ time_point)
$emmeans
time_point emmean SE df asymp.LCL asymp.UCL
Acclimation 4.05 0.0938 Inf 3.87 4.24
FMT1 3.84 0.0949 Inf 3.66 4.03
FMT2 4.21 0.0895 Inf 4.03 4.38
Results are averaged over the levels of: type
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Acclimation - FMT1 0.209 0.123 Inf 1.700 0.2050
Acclimation - FMT2 -0.156 0.121 Inf -1.286 0.4028
FMT1 - FMT2 -0.365 0.122 Inf -2.997 0.0077
Results are averaged over the levels of: type
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 3 estimates
Neutral
Modelq1n <- nlme::lme(neutral ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Modelq1n)
Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
563.2664 587.9998 -270.6332
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 4.551076 9.160157
Fixed effects: neutral ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 21.343055 3.603126 46 5.923483 0.0000
typeHot_control 23.258307 4.960549 24 4.688656 0.0001
typeTreatment -3.944072 4.960549 24 -0.795088 0.4344
time_pointFMT1 -3.983722 4.472618 46 -0.890691 0.3777
time_pointFMT2 2.374489 4.472618 46 0.530895 0.5980
typeHot_control:time_pointFMT1 -14.877819 6.216964 46 -2.393101 0.0208
typeTreatment:time_pointFMT1 11.253000 6.325237 46 1.779064 0.0818
typeHot_control:time_pointFMT2 -14.775497 6.216964 46 -2.376642 0.0217
typeTreatment:time_pointFMT2 20.108022 6.216964 46 3.234380 0.0023
Correlation:
(Intr) typHt_ typTrt t_FMT1 t_FMT2 tH_:_FMT1 tT:_FMT1 tH_:_FMT2
typeHot_control -0.726
typeTreatment -0.726 0.528
time_pointFMT1 -0.663 0.481 0.481
time_pointFMT2 -0.663 0.481 0.481 0.534
typeHot_control:time_pointFMT1 0.477 -0.649 -0.346 -0.719 -0.384
typeTreatment:time_pointFMT1 0.469 -0.340 -0.638 -0.707 -0.378 0.509
typeHot_control:time_pointFMT2 0.477 -0.649 -0.346 -0.384 -0.719 0.518 0.272
typeTreatment:time_pointFMT2 0.477 -0.346 -0.649 -0.384 -0.719 0.276 0.509 0.518
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.1324340 -0.7088533 0.0456110 0.6569569 1.9681295
Number of Observations: 79
Number of Groups: 27
emmeans::emmeans(Modelq1n, pairwise ~ type)
$emmeans
type emmean SE df lower.CL upper.CL
Control 20.8 2.36 26 16.0 25.7
Hot_control 34.2 2.33 24 29.4 39.0
Treatment 27.3 2.36 24 22.4 32.2
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Control - Hot_control -13.37 3.31 24 -4.038 0.0013
Control - Treatment -6.51 3.33 24 -1.952 0.1461
Hot_control - Treatment 6.86 3.31 24 2.073 0.1170
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates
emmeans::emmeans(Modelq1n, pairwise ~ time_point)
$emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 27.8 2.01 24 23.6 31.9
FMT1 22.6 2.01 24 18.4 26.7
FMT2 31.9 1.97 24 27.9 36.0
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - FMT1 5.19 2.55 46 2.034 0.1156
Acclimation - FMT2 -4.15 2.52 46 -1.646 0.2372
FMT1 - FMT2 -9.34 2.52 46 -3.703 0.0016
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates
Phylogenetic
Model_phylo <- nlme::lme(phylogenetic ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Model_phylo)
Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
269.7976 294.5311 -123.8988
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.4986226 1.145917
Fixed effects: phylogenetic ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 5.256251 0.4407611 46 11.925398 0.0000
typeHot_control 1.266332 0.6064636 24 2.088060 0.0476
typeTreatment 0.283689 0.6064636 24 0.467776 0.6442
time_pointFMT1 -0.837008 0.5590601 46 -1.497169 0.1412
time_pointFMT2 -0.481515 0.5590601 46 -0.861293 0.3935
typeHot_control:time_pointFMT1 -1.410952 0.7774021 46 -1.814958 0.0761
typeTreatment:time_pointFMT1 -0.545877 0.7906304 46 -0.690432 0.4934
typeHot_control:time_pointFMT2 -0.585100 0.7774021 46 -0.752635 0.4555
typeTreatment:time_pointFMT2 0.056099 0.7774021 46 0.072162 0.9428
Correlation:
(Intr) typHt_ typTrt t_FMT1 t_FMT2 tH_:_FMT1 tT:_FMT1 tH_:_FMT2
typeHot_control -0.727
typeTreatment -0.727 0.528
time_pointFMT1 -0.676 0.492 0.492
time_pointFMT2 -0.676 0.492 0.492 0.533
typeHot_control:time_pointFMT1 0.486 -0.663 -0.353 -0.719 -0.383
typeTreatment:time_pointFMT1 0.478 -0.348 -0.652 -0.707 -0.377 0.509
typeHot_control:time_pointFMT2 0.486 -0.663 -0.353 -0.383 -0.719 0.517 0.271
typeTreatment:time_pointFMT2 0.486 -0.353 -0.663 -0.383 -0.719 0.276 0.509 0.517
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.81722823 -0.52018811 -0.06481484 0.55835817 2.23628404
Number of Observations: 79
Number of Groups: 27
emmeans::emmeans(Model_phylo, pairwise ~ type)
$emmeans
type emmean SE df lower.CL upper.CL
Control 4.82 0.280 26 4.24 5.39
Hot_control 5.42 0.276 24 4.85 5.99
Treatment 4.94 0.280 24 4.36 5.52
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Control - Hot_control -0.601 0.393 24 -1.527 0.2963
Control - Treatment -0.120 0.396 24 -0.304 0.9505
Hot_control - Treatment 0.481 0.393 24 1.221 0.4525
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates
emmeans::emmeans(Model_phylo, pairwise ~ time_point)
$emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 5.77 0.245 24 5.27 6.28
FMT1 4.28 0.245 24 3.78 4.79
FMT2 5.12 0.241 24 4.62 5.61
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - FMT1 1.489 0.319 46 4.666 0.0001
Acclimation - FMT2 0.658 0.316 46 2.085 0.1042
FMT1 - FMT2 -0.831 0.316 46 -2.635 0.0302
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates
Beta diversity
Number of samples used
<<<<<<< HEAD
samples_to_keep_post7 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation"|time_point=="FMT2") %>%
select(Tube_code) %>%
pull()
subset_meta_post7 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation"|time_point=="FMT2")
subset_meta_post7$time_point<-as.factor(subset_meta_post7$time_point)
subset_meta_post7$type<-as.factor(subset_meta_post7$type)
length(samples_to_keep_post7)
[1] 79
Richness
richness_post7 <- as.matrix(beta_q0n$S)
richness_post7 <- as.dist(richness_post7[rownames(richness_post7) %in% samples_to_keep_post7,
colnames(richness_post7) %in% samples_to_keep_post7])
betadisper(richness_post7, subset_meta_post7$type) %>% permutest(., pairwise = TRUE)
=======
samples_to_keep_post3 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Hot_control") %>%
select(Tube_code) %>%
pull()
subset_meta_post3 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Hot_control")
subset_meta_post3$type_time<-interaction(subset_meta_post3$type, subset_meta_post3$time_point)
length(samples_to_keep_post3)
[1] 34
Richness
richness_post3 <- as.matrix(beta_q0n$S)
richness_post3 <- as.dist(richness_post3[rownames(richness_post3) %in% samples_to_keep_post3,
colnames(richness_post3) %in% samples_to_keep_post3])
betadisper(richness_post3, subset_meta_post3$type) %>% permutest(., pairwise = TRUE)
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.06454 0.032268 4.9753 999 0.01 **
Residuals 76 0.49290 0.006485
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
<<<<<<< HEAD
Control Hot_control Treatment
Control 0.0160000 0.697
Hot_control 0.0119513 0.002
Treatment 0.6817089 0.0035106
adonis2(richness_post7 ~ type*time_point,
data = subset_meta_post7 %>% arrange(match(Tube_code,labels(richness_post7))),
permutations = 999,
strata = subset_meta_post7 %>% arrange(match(Tube_code,labels(richness_post7))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()
=======
Control Treatment
Control 0.484
Treatment 0.49714
adonis2(richness_post3 ~ type*time_point,
data = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))),
permutations = 999,
strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
=======
.table td.tinytable_css_b9npyawa3sc60nzpd9sg, .table th.tinytable_css_b9npyawa3sc60nzpd9sg { border-bottom: solid #d3d8dc 0.1em; }
.table td.tinytable_css_375734g7wmtqf6o2enbc, .table th.tinytable_css_375734g7wmtqf6o2enbc { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| type |
2 |
3.120806 |
0.11961754 |
5.546155 |
0.001 |
| time_point |
2 |
1.703776 |
0.06530413 |
3.027874 |
0.001 |
| type:time_point |
4 |
1.570884 |
0.06021051 |
1.395852 |
0.001 |
| Residual |
70 |
19.694403 |
0.75486782 |
NA |
NA |
| Total |
78 |
26.089870 |
1.00000000 |
NA |
NA |
<<<<<<< HEAD
subset_meta_post7_arrange <- column_to_rownames(subset_meta_post7, "Tube_code")
subset_meta_post7_arrange<-subset_meta_post7_arrange[labels(richness_post7),]
pairwise <- pairwise.adonis(richness_post7, subset_meta_post7_arrange$type, perm=999)
pairwise%>%
tt()
subset_meta_post3_arrange <- column_to_rownames(subset_meta_post3, "Tube_code")
subset_meta_post3_arrange<-subset_meta_post3_arrange[labels(richness_post3),]
pairwise <- pairwise.adonis(richness_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
tt()
=======
.table td.tinytable_css_tbs2p2aupt3txp1lbmlr, .table th.tinytable_css_tbs2p2aupt3txp1lbmlr { border-bottom: solid #d3d8dc 0.1em; }
.table td.tinytable_css_a6pij8s035i9eixvhixm, .table th.tinytable_css_a6pij8s035i9eixvhixm { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Control vs Treatment |
1 |
<<<<<<< HEAD
0.7820564 |
2.401526 |
0.04582931 |
0.001 |
0.003 |
* |
=======
0.3657243 |
1.123239 |
0.06966584 |
0.248 |
1.000 |
|
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| Control vs Hot_control |
1 |
2.6071063 |
9.057490 |
0.15081366 |
0.001 |
0.003 |
* |
| Treatment vs Hot_control |
1 |
1.2773606 |
4.350043 |
0.07859150 |
0.001 |
0.003 |
* |
pairwise <- pairwise.adonis(richness_post7, subset_meta_post7_arrange$time_point, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Acclimation vs FMT1 |
1 |
1.1985353 |
3.576695 |
0.06675842 |
0.001 |
0.003 |
* |
| Acclimation vs FMT2 |
1 |
0.7814316 |
2.468536 |
0.04616802 |
0.003 |
0.009 |
* |
| FMT1 vs FMT2 |
1 |
<<<<<<< HEAD
0.5620083 |
1.802651 |
0.03413940 |
0.011 |
0.033 |
. |
=======
0.5615418 |
1.729004 |
0.10335366 |
0.017 |
0.102 |
|
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
Neutral
<<<<<<< HEAD
neutral_post7 <- as.matrix(beta_q1n$S)
neutral_post7 <- as.dist(neutral_post7[rownames(neutral_post7) %in% samples_to_keep_post7,
colnames(neutral_post7) %in% samples_to_keep_post7])
betadisper(neutral_post7, subset_meta_post7$type) %>% permutest(., pairwise = TRUE)
neutral_post3 <- as.matrix(beta_q1n$S)
neutral_post3 <- as.dist(neutral_post3[rownames(neutral_post3) %in% samples_to_keep_post3,
colnames(neutral_post3) %in% samples_to_keep_post3])
betadisper(neutral_post3, subset_meta_post3$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
<<<<<<< HEAD
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.06738 0.033689 3.8309 999 0.024 *
Residuals 76 0.66833 0.008794
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.2060000 0.173
Hot_control 0.1985046 0.005
Treatment 0.1611083 0.0069718
adonis2(neutral_post7 ~ type*time_point,
data = subset_meta_post7 %>% arrange(match(Tube_code,labels(neutral_post7))),
permutations = 999,
strata = subset_meta_post7 %>% arrange(match(Tube_code,labels(neutral_post7))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()
adonis2(neutral_post3 ~ type,
data = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))),
permutations = 999,
strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()
=======
.table td.tinytable_css_z8wsgq48zbrihjwlnh7p, .table th.tinytable_css_z8wsgq48zbrihjwlnh7p { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
.table td.tinytable_css_dj3n5nfx6r2y9emr2fzn, .table th.tinytable_css_dj3n5nfx6r2y9emr2fzn { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
<<<<<<< HEAD
| type |
2 |
3.705559 |
0.14953527 |
7.609355 |
0.001 |
| time_point |
2 |
2.217656 |
0.08949199 |
4.553951 |
0.001 |
| type:time_point |
4 |
1.813192 |
0.07317011 |
1.861693 |
0.001 |
=======
Model |
1 |
0.4361022 |
0.04050485 |
1.350872 |
1 |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| Residual |
70 |
17.044094 |
0.68780262 |
NA |
NA |
| Total |
78 |
24.780501 |
1.00000000 |
NA |
NA |
<<<<<<< HEAD
subset_meta_post7_arrange <- column_to_rownames(subset_meta_post7, "Tube_code")
subset_meta_post7_arrange<-subset_meta_post7_arrange[labels(neutral_post7),]
pairwise <- pairwise.adonis(neutral_post7, subset_meta_post7$type, perm=999)
pairwise%>%
tt()
subset_meta_post3_arrange <- column_to_rownames(subset_meta_post3, "Tube_code")
subset_meta_post3_arrange<-subset_meta_post3_arrange[labels(neutral_post3),]
pairwise <- pairwise.adonis(neutral_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
tt()
=======
.table td.tinytable_css_iagsoasjtf3ipwhzulgy, .table th.tinytable_css_iagsoasjtf3ipwhzulgy { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
.table td.tinytable_css_a18ysd90ugpwhxe6heg0, .table th.tinytable_css_a18ysd90ugpwhxe6heg0 { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Control vs Treatment |
1 |
<<<<<<< HEAD
1.027801 |
3.463893 |
0.06478939 |
=======
0.2759858 |
0.9928976 |
0.06208366 |
0.483 |
1.000 |
|
| Control.Acclimation vs Control.FMT1 |
1 |
0.8153894 |
3.0970603 |
0.17113610 |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
0.001 |
0.003 |
* |
| Control vs Hot_control |
1 |
<<<<<<< HEAD
3.120025 |
12.118158 |
0.19199163 |
0.001 |
0.003 |
* |
=======
1.1809241 |
4.4856470 |
0.24265567 |
0.002 |
0.012 |
. |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| Treatment vs Hot_control |
1 |
1.394947 |
5.015972 |
0.08954539 |
0.001 |
0.003 |
* |
<<<<<<< HEAD
=======
| Treatment.Acclimation vs Treatment.FMT1 |
1 |
1.3127773 |
4.6298256 |
0.23585668 |
0.001 |
0.006 |
* |
| Control.FMT1 vs Treatment.FMT1 |
1 |
0.6051778 |
2.2508491 |
0.13047758 |
0.015 |
0.090 |
|
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
<<<<<<< HEAD
pairwise <- pairwise.adonis(neutral_post7, subset_meta_post7_arrange$time_point, perm=999)
pairwise%>%
tt()
Phylogenetic
phylo_post3 <- as.matrix(beta_q1p$S)
phylo_post3 <- as.dist(phylo_post3[rownames(phylo_post3) %in% samples_to_keep_post3,
colnames(phylo_post3) %in% samples_to_keep_post3])
betadisper(phylo_post3, subset_meta_post3$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04697 0.046974 2.8076 999 0.101
Residuals 32 0.53540 0.016731
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Treatment
Control 0.096
Treatment 0.10357
adonis2(phylo_post3 ~ type,
data = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))),
permutations = 999,
strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()
=======
.table td.tinytable_css_khvi88bgx62bf4888wka, .table th.tinytable_css_khvi88bgx62bf4888wka { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
.table td.tinytable_css_5y8xbfi1thjnoay3akf7, .table th.tinytable_css_5y8xbfi1thjnoay3akf7 { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Acclimation vs FMT1 |
1 |
<<<<<<< HEAD
1.6914939 |
5.559526 |
0.10006431 |
0.001 |
0.003 |
* |
=======
0.03710337 |
0.01720754 |
0.5602822 |
1 |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| Acclimation vs FMT2 |
1 |
0.9608387 |
3.212301 |
0.05925409 |
0.001 |
0.003 |
* |
| FMT1 vs FMT2 |
1 |
0.6282036 |
2.174336 |
0.04089070 |
0.010 |
0.030 |
. |
<<<<<<< HEAD
Phylogenetic
phylo_post7 <- as.matrix(beta_q1p$S)
phylo_post7 <- as.dist(phylo_post7[rownames(phylo_post7) %in% samples_to_keep_post7,
colnames(phylo_post7) %in% samples_to_keep_post7])
betadisper(phylo_post7, subset_meta_post7$type) %>% permutest(., pairwise = TRUE)
NMDS
beta_richness_nmds_post3 <- richness_post3 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_post3, by = join_by(sample == Tube_code)) %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
beta_neutral_nmds_post3 <- neutral_post3 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_post3, by = join_by(sample == Tube_code))%>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
beta_phylogenetic_nmds_post3 <- phylo_post3 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_post3, by = join_by(sample == Tube_code))%>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
ggarrange(beta_richness_nmds_post3, beta_neutral_nmds_post3, beta_phylogenetic_nmds_post3, ncol=2, nrow=2, common.legend = TRUE, legend="right")
Functional differences
CC from acclimation to FMT1
element_gift %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Control") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 11 × 6
trait p_value p_adjust Domain Function Element
<chr> <dbl> <dbl> <chr> <chr> <chr>
1 B0204 0.00370 0.0478 Biosynthesis Amino acid biosynthesis_Serine Serine
2 B0205 0.000165 0.0117 Biosynthesis Amino acid biosynthesis_Threonine Threonine
3 B0703 0.00370 0.0478 Biosynthesis Vitamin biosynthesis_Niacin (B3) Niacin (B3)
4 B0805 0.00156 0.0370 Biosynthesis Aromatic compound biosynthesis_Indole-3-acetate Indole-3-acetate
5 D0102 0.000329 0.0156 Degradation Lipid degradation_Fatty acid Fatty acid
6 D0212 0.000987 0.0350 Degradation Polysaccharide degradation_Arabinan Arabinan
7 D0304 0.00247 0.0438 Degradation Sugar degradation_D-Arabinose D-Arabinose
8 D0306 0.00156 0.0370 Degradation Sugar degradation_D-Xylose D-Xylose
9 D0308 0.00247 0.0438 Degradation Sugar degradation_L-Rhamnose L-Rhamnose
10 D0601 0.00370 0.0478 Degradation Nitrogen compound degradation_Nitrate Nitrate
11 D0708 0.000165 0.0117 Degradation Alcohol degradation_Phytol Phytol
CI from acclimation to FMT1
element_gift %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Treatment") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 2 × 6
trait p_value p_adjust Domain Function Element
<chr> <dbl> <dbl> <chr> <chr> <chr>
1 B0710 0.0000823 0.0119 Biosynthesis Vitamin biosynthesis_Phylloquinone (K1) Phylloquinone (K1)
2 D0102 0.000576 0.0418 Degradation Lipid degradation_Fatty acid Fatty acid
CI vs WC
Alpha
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Control") %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
geom_boxplot(outlier.shape = NA) +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5) +
facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
values=c("#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
values=c("#d57d2c50","#76b18350")) +
theme_minimal() +
theme(axis.text.x=element_text(size=8))+
labs(x = "Time Point", y = "value")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Control")
Richness
Modelq0GLMMNB <- MASS::glm.nb(richness ~ type*time_point, data = alpha_div_meta,trace=TRUE)
summary(Modelq0GLMMNB)
Call:
MASS::glm.nb(formula = richness ~ type * time_point, data = alpha_div_meta,
trace = TRUE, init.theta = 5.042866344, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.4697 0.1527 29.279 < 2e-16 ***
typeTreatment -0.6557 0.2186 -2.999 0.00271 **
time_pointFMT1 -0.4209 0.2174 -1.936 0.05292 .
typeTreatment:time_pointFMT1 0.5484 0.3146 1.743 0.08131 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(5.0429) family taken to be 1)
Null deviance: 48.066 on 34 degrees of freedom
Residual deviance: 37.426 on 31 degrees of freedom
AIC: 336.8
Number of Fisher Scoring iterations: 1
Theta: 5.04
Std. Err.: 1.34
2 x log-likelihood: -326.796
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ type)
$emmeans
type emmean SE df asymp.LCL asymp.UCL
Hot_control 4.26 0.109 Inf 4.05 4.47
Treatment 3.88 0.114 Inf 3.65 4.10
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Hot_control - Treatment 0.381 0.157 Inf 2.425 0.0153
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Neutral
Modelq1n <- lm(formula = neutral ~ type*time_point, data = alpha_div_meta)
summary(Modelq1n)
Call:
lm(formula = neutral ~ type * time_point, data = alpha_div_meta)
Residuals:
Min 1Q Median 3Q Max
-20.6988 -6.3363 -0.4774 7.6008 15.8281
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 44.601 3.324 13.418 1.88e-14 ***
typeTreatment -27.202 4.701 -5.786 2.26e-06 ***
time_pointFMT1 -18.862 4.701 -4.012 0.000353 ***
typeTreatment:time_pointFMT1 26.326 6.751 3.899 0.000483 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.972 on 31 degrees of freedom
Multiple R-squared: 0.5397, Adjusted R-squared: 0.4951
F-statistic: 12.12 on 3 and 31 DF, p-value: 2.052e-05
Phylogenetic
Model_phylo <- lm(formula = phylogenetic ~ type*time_point, data = alpha_div_meta)
summary(Model_phylo)
Call:
lm(formula = phylogenetic ~ type * time_point, data = alpha_div_meta)
Residuals:
Min 1Q Median 3Q Max
-3.07987 -0.70068 -0.07723 0.54734 2.84255
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.5226 0.4274 15.260 5.8e-16 ***
typeTreatment -0.9826 0.6045 -1.626 0.114169
time_pointFMT1 -2.2480 0.6045 -3.719 0.000793 ***
typeTreatment:time_pointFMT1 0.8981 0.8681 1.035 0.308891
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.282 on 31 degrees of freedom
Multiple R-squared: 0.39, Adjusted R-squared: 0.331
F-statistic: 6.606 on 3 and 31 DF, p-value: 0.001398
Beta diversity
Number of samples used
samples_to_keep_post4 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Control") %>%
select(Tube_code) %>%
pull()
subset_meta_post4 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Control")
subset_meta_post4$type_time<-interaction(subset_meta_post4$type, subset_meta_post4$time_point)
length(samples_to_keep_post4)
[1] 35
Richness
richness_post4 <- as.matrix(beta_q0n$S)
richness_post4 <- as.dist(richness_post4[rownames(richness_post4) %in% samples_to_keep_post4,
colnames(richness_post4) %in% samples_to_keep_post4])
betadisper(richness_post4, subset_meta_post4$type) %>% permutest(., pairwise = TRUE)
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.05438 0.027192 2.1499 999 0.126
Residuals 76 0.96123 0.012648
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
<<<<<<< HEAD
Control Hot_control Treatment
Control 0.718000 0.160
Hot_control 0.697512 0.084
Treatment 0.155309 0.083221
adonis2(phylo_post7 ~ type*time_point,
data = subset_meta_post7 %>% arrange(match(Tube_code,labels(phylo_post7))),
permutations = 999,
strata = subset_meta_post7 %>% arrange(match(Tube_code,labels(phylo_post7))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()
=======
Hot_control Treatment
Hot_control 0.001
Treatment 0.00010134
adonis2(richness_post4 ~ type*time_point,
data = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))),
permutations = 999,
strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
=======
.table td.tinytable_css_zyvbojr9jfe6h7isy3qs, .table th.tinytable_css_zyvbojr9jfe6h7isy3qs { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
.table td.tinytable_css_q3vlmkcbossuw9eaeyjx, .table th.tinytable_css_q3vlmkcbossuw9eaeyjx { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| type |
2 |
0.2941494 |
0.07500331 |
3.592254 |
0.001 |
| time_point |
2 |
0.5528857 |
0.14097688 |
6.752032 |
0.001 |
| type:time_point |
4 |
0.2088312 |
0.05324856 |
1.275159 |
0.152 |
| Residual |
70 |
2.8659522 |
0.73077126 |
NA |
NA |
| Total |
78 |
3.9218184 |
1.00000000 |
NA |
NA |
<<<<<<< HEAD
subset_meta_post7_arrange <- column_to_rownames(subset_meta_post7, "Tube_code")
subset_meta_post7_arrange<-subset_meta_post7_arrange[labels(phylo_post7),]
pairwise <- pairwise.adonis(phylo_post7, subset_meta_post7_arrange$type, perm=999)
pairwise%>%
tt()
subset_meta_post4_arrange <- column_to_rownames(subset_meta_post4, "Tube_code")
subset_meta_post4_arrange<-subset_meta_post4_arrange[labels(richness_post4),]
pairwise <- pairwise.adonis(richness_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
tt()
=======
.table td.tinytable_css_xdkcihx41r32vi6u5eip, .table th.tinytable_css_xdkcihx41r32vi6u5eip { border-bottom: solid #d3d8dc 0.1em; }
.table td.tinytable_css_4l2ubp19yc5gisw70e31, .table th.tinytable_css_4l2ubp19yc5gisw70e31 { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
<<<<<<< HEAD
| Control vs Treatment |
=======
Treatment.Acclimation vs Hot_control.Acclimation |
1 |
1.3606630 |
5.087152 |
0.2412441 |
0.002 |
0.012 |
. |
| Treatment.Acclimation vs Treatment.FMT1 |
1 |
0.9551308 |
2.926054 |
0.1632291 |
0.001 |
0.006 |
* |
| Treatment.Acclimation vs Hot_control.FMT1 |
1 |
1.2263345 |
4.039487 |
0.2015764 |
0.001 |
0.006 |
* |
| Hot_control.Acclimation vs Treatment.FMT1 |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
1 |
0.03324557 |
0.6023122 |
0.01190286 |
0.649 |
1.000 |
|
| Control vs Hot_control |
1 |
0.22899071 |
6.3999219 |
0.11149705 |
0.001 |
0.003 |
* |
| Treatment vs Hot_control |
1 |
<<<<<<< HEAD
0.17684005 |
3.3769360 |
0.06210236 |
0.014 |
0.042 |
. |
=======
0.3734921 |
1.268929 |
0.0779971 |
0.110 |
0.660 |
|
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
<<<<<<< HEAD
pairwise <- pairwise.adonis(phylo_post7, subset_meta_post7_arrange$time_point, perm=999)
pairwise%>%
tt()
Neutral
neutral_post4 <- as.matrix(beta_q1n$S)
neutral_post4 <- as.dist(neutral_post4[rownames(neutral_post4) %in% samples_to_keep_post4,
colnames(neutral_post4) %in% samples_to_keep_post4])
betadisper(neutral_post4, subset_meta_post4$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.087335 0.087335 15.924 999 0.003 **
Residuals 33 0.180987 0.005484
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
Hot_control 0.003
Treatment 0.00034565
adonis2(neutral_post4 ~ type,
data = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))),
permutations = 999,
strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()
=======
.table td.tinytable_css_xxa6ri3auop2jfg0ti55, .table th.tinytable_css_xxa6ri3auop2jfg0ti55 { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
.table td.tinytable_css_vwsupixy2hc38u228396, .table th.tinytable_css_vwsupixy2hc38u228396 { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Acclimation vs FMT1 |
1 |
<<<<<<< HEAD
0.3794332 |
7.332180 |
0.12788943 |
0.001 |
0.003 |
* |
=======
1.332543 |
0.1214599 |
4.562313 |
1 |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| Acclimation vs FMT2 |
1 |
0.3351380 |
7.390585 |
0.12657152 |
0.001 |
0.003 |
* |
| FMT1 vs FMT2 |
1 |
0.1178756 |
3.274037 |
0.06032418 |
0.005 |
0.015 |
. |
<<<<<<< HEAD
dbRDA
#Richness
cca_ord <- capscale(formula = richness_post7 ~ subset_meta_post7$time_point* subset_meta_post7$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post7, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post7$time_pointFMT1" = "FMT1",
"subset_meta_post7$time_pointFMT2"="FMT2",
"subset_meta_post7$typeHot_control"="Warm-control",
"subset_meta_post7$typeTreatment" = "Cold-intervention",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeHot_control"="FMT1 Warm-control",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeHot_control"="FMT2 Warm-control",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeTreatment" = "FMT1 Cold-intervention",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeTreatment" = "FMT2 Cold-intervention")
beta_richness_nmds_post7 <-CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-control","Warm-control", "Cold-intervention"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-control","Warm-control", "Cold-intervention"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Neutral
cca_ord <- capscale(formula = neutral_post7 ~ subset_meta_post7$time_point* subset_meta_post7$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post7, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post7$time_pointFMT1" = "FMT1",
"subset_meta_post7$time_pointFMT2"="FMT2",
"subset_meta_post7$typeHot_control"="Warm-control",
"subset_meta_post7$typeTreatment" = "Cold-intervention",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeHot_control"="FMT1 Warm-control",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeHot_control"="FMT2 Warm-control",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeTreatment" = "FMT1 Cold-intervention",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeTreatment" = "FMT2 Cold-intervention")
beta_neutral_nmds_post7 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-control","Warm-control", "Cold-intervention"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-control","Warm-control", "Cold-intervention"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Phylogenetic
cca_ord <- capscale(formula = phylo_post7 ~ subset_meta_post7$time_point* subset_meta_post7$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post7, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post7$time_pointFMT1" = "FMT1",
"subset_meta_post7$time_pointFMT2"="FMT2",
"subset_meta_post7$typeHot_control"="Warm-control",
"subset_meta_post7$typeTreatment" = "Cold-intervention",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeHot_control"="FMT1 Warm-control",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeHot_control"="FMT2 Warm-control",
"subset_meta_post7$time_pointFMT1:subset_meta_post7$typeTreatment" = "FMT1 Cold-intervention",
"subset_meta_post7$time_pointFMT2:subset_meta_post7$typeTreatment" = "FMT2 Cold-intervention")
beta_phylogenetic_nmds_post7 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-control","Warm-control", "Cold-intervention"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-control","Warm-control", "Cold-intervention"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
ggarrange(beta_richness_nmds_post7, beta_neutral_nmds_post7, beta_phylogenetic_nmds_post7, ncol=3, nrow=1, common.legend = TRUE, legend="right")
What is the effect of FMT and temperature on the GM after 1 week compared to acclimation?
CI vs CC
Alpha diversity
label_map <- c(
"Control" = "Cold-control",
"Treatment" = "Cold-intervention",
"richness" = "Species Richness",
"neutral" = "Neutral Diversity",
"phylogenetic" = "Phylogenetic Diversity"
)
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(type %in% c("Control","Treatment") & time_point %in% c("FMT1","Acclimation")) %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5, show.legend = FALSE) +
facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
# facet_nested(.~metric+type, labeller = labeller(metric = label_map, type = label_map))+
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
values=c("#4477AA", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control", "Treatment"),
values=c("#4477AA50","#76b18350")) +
theme_minimal() +
theme(axis.text.x=element_text(size=10))+
labs(x = "Time Point", y = "value")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Hot_control")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(2.5327) ( log )
Formula: richness ~ type * time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
329.2 338.4 -158.6 317.2 28
Scaled residuals:
Min 1Q Median 3Q Max
-1.4046 -0.7191 0.1331 0.6747 1.4330
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 2.062e-11 4.541e-06
Number of obs: 34, groups: individual, 18
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.9512 0.2275 17.368 <2e-16 ***
typeTreatment -0.1372 0.3132 -0.438 0.661
time_pointFMT1 -0.3225 0.3140 -1.027 0.304
typeTreatment:time_pointFMT1 0.4500 0.4435 1.015 0.310
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) typTrt t_FMT1
typeTretmnt -0.726
tim_pntFMT1 -0.725 0.526
typTr:_FMT1 0.513 -0.706 -0.708
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ type)
$emmeans
type emmean SE df asymp.LCL asymp.UCL
Control 3.79 0.157 Inf 3.48 4.10
Treatment 3.88 0.157 Inf 3.57 4.18
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Control - Treatment -0.0878 0.222 Inf -0.396 0.6921
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ time_point)
$emmeans
time_point emmean SE df asymp.LCL asymp.UCL
Acclimation 3.88 0.157 Inf 3.58 4.19
FMT1 3.79 0.157 Inf 3.48 4.09
Results are averaged over the levels of: type
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Acclimation - FMT1 0.0975 0.222 Inf 0.440 0.6603
Results are averaged over the levels of: type
Results are given on the log (not the response) scale.
Neutral
Modelq1n <- lme(neutral ~ type*time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Modelq1n)
Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
244.9928 253.4 -116.4964
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 6.116057 8.521967
Fixed effects: neutral ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 21.027978 3.684717 16 5.706810 0.0000
typeTreatment -3.628995 5.079636 16 -0.714420 0.4853
time_pointFMT1 -3.668645 4.182131 14 -0.877219 0.3952
typeTreatment:time_pointFMT1 10.917012 5.914427 14 1.845828 0.0862
Correlation:
(Intr) typTrt t_FMT1
typeTreatment -0.725
time_pointFMT1 -0.611 0.443
typeTreatment:time_pointFMT1 0.432 -0.582 -0.707
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.79785802 -0.67295800 0.08387402 0.71606547 1.30935708
Number of Observations: 34
Number of Groups: 18
emmeans::emmeans(Modelq1n, pairwise ~ type)
$emmeans
type emmean SE df lower.CL upper.CL
Control 19.2 2.92 14 12.9 25.5
Treatment 21.0 2.92 14 14.8 27.3
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Control - Treatment -1.83 4.13 14 -0.443 0.6646
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
emmeans::emmeans(Modelq1n, pairwise ~ time_point)
$emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 19.2 2.54 16 13.8 24.6
FMT1 21.0 2.54 14 15.6 26.5
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - FMT1 -1.79 2.96 14 -0.605 0.5547
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Phylogenetic
Model_phylo <- lme(phylogenetic ~ type*time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)
Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
125.282 133.6892 -56.64101
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 7.854299e-05 1.386163
Fixed effects: phylogenetic ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 5.227341 0.4900827 16 10.666243 0.0000
typeTreatment 0.312600 0.6735542 16 0.464105 0.6488
time_pointFMT1 -0.808097 0.6735542 14 -1.199751 0.2501
typeTreatment:time_pointFMT1 -0.541741 0.9525495 14 -0.568728 0.5786
Correlation:
(Intr) typTrt t_FMT1
typeTreatment -0.728
time_pointFMT1 -0.728 0.529
typeTreatment:time_pointFMT1 0.514 -0.707 -0.707
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.630942795 -0.368275262 0.004976243 0.364227875 2.050663052
Number of Observations: 34
Number of Groups: 18
emmeans::emmeans(Model_phylo, pairwise ~ type)
$emmeans
type emmean SE df lower.CL upper.CL
Control 4.82 0.337 14 4.10 5.55
Treatment 4.87 0.337 14 4.14 5.59
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Control - Treatment -0.0417 0.476 14 -0.088 0.9314
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
emmeans::emmeans(Model_phylo, pairwise ~ time_point)
$emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 5.38 0.337 16 4.67 6.10
FMT1 4.30 0.337 14 3.58 5.03
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - FMT1 1.08 0.476 14 2.265 0.0399
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Beta diversity
Number of samples used
samples_to_keep_post3 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Hot_control") %>%
select(Tube_code) %>%
pull()
subset_meta_post3 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Hot_control")
subset_meta_post3$time_point<-as.factor(subset_meta_post3$time_point)
subset_meta_post3$type<-as.factor(subset_meta_post3$type)
length(samples_to_keep_post3)
[1] 34
Richness
richness_post3 <- as.matrix(beta_q0n$S)
richness_post3 <- as.dist(richness_post3[rownames(richness_post3) %in% samples_to_keep_post3,
colnames(richness_post3) %in% samples_to_keep_post3])
betadisper(richness_post3, subset_meta_post3$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.002424 0.0024241 0.4717 999 0.482
Residuals 32 0.164445 0.0051389
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Treatment
Control 0.47
Treatment 0.49714
adonis2(richness_post3 ~ type*time_point,
data = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))),
permutations = 999,
strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()
=======
subset_meta_post4_arrange <- column_to_rownames(subset_meta_post4, "Tube_code")
subset_meta_post4_arrange<-subset_meta_post4_arrange[labels(neutral_post4),]
pairwise <- pairwise.adonis(neutral_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
=======
.table td.tinytable_css_xequglxj0vva2o6kdxtz, .table th.tinytable_css_xequglxj0vva2o6kdxtz { border-bottom: solid #d3d8dc 0.1em; }
.table td.tinytable_css_hxyn4amlq8uzaste8w29, .table th.tinytable_css_hxyn4amlq8uzaste8w29 { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| type |
1 |
0.4269704 |
0.03612811 |
1.312996 |
0.001 |
| time_point |
1 |
1.1427754 |
0.09669595 |
3.514201 |
0.001 |
| type:time_point |
1 |
0.4928541 |
0.04170286 |
1.515598 |
0.050 |
| Residual |
30 |
9.7556341 |
0.82547308 |
NA |
NA |
<<<<<<< HEAD
| Total |
33 |
11.8182341 |
1.00000000 |
NA |
NA |
=======
Treatment.FMT1 vs Hot_control.FMT1 |
1 |
0.4150076 |
1.637268 |
0.09840968 |
0.062 |
0.372 |
|
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
<<<<<<< HEAD
subset_meta_post3_arrange <- column_to_rownames(subset_meta_post3, "Tube_code")
subset_meta_post3_arrange<-subset_meta_post3_arrange[labels(richness_post3),]
pairwise <- pairwise.adonis(richness_post3, subset_meta_post3_arrange$type, perm=999)
pairwise%>%
tt()
Phylogenetic
phylo_post4 <- as.matrix(beta_q1p$S)
phylo_post4 <- as.dist(phylo_post4[rownames(phylo_post4) %in% samples_to_keep_post4,
colnames(phylo_post4) %in% samples_to_keep_post4])
betadisper(phylo_post4, subset_meta_post4$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.07659 0.076588 5.0159 999 0.026 *
Residuals 33 0.50388 0.015269
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
Hot_control 0.033
Treatment 0.031971
adonis2(phylo_post4 ~ type,
data = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))),
permutations = 999,
strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()
=======
.table td.tinytable_css_l5js7rbeqaxhumvczvgh, .table th.tinytable_css_l5js7rbeqaxhumvczvgh { border-bottom: solid #d3d8dc 0.1em; }
.table td.tinytable_css_k4fj9c6v68b3z2lu952e, .table th.tinytable_css_k4fj9c6v68b3z2lu952e { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Control vs Treatment |
1 |
<<<<<<< HEAD
0.4269704 |
1.199433 |
0.03612811 |
0.166 |
0.166 |
|
=======
0.2070779 |
0.09373244 |
3.413088 |
1 |
| Residual |
33 |
2.0021670 |
0.90626756 |
NA |
NA |
| Total |
34 |
2.2092449 |
1.00000000 |
NA |
NA |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
<<<<<<< HEAD
pairwise <- pairwise.adonis(richness_post3, subset_meta_post3_arrange$time_point, perm=999)
pairwise%>%
tt()
NMDS
beta_richness_nmds_post4 <- richness_post4 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_post4, by = join_by(sample == Tube_code)) %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
beta_neutral_nmds_post4 <- neutral_post4 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_post4, by = join_by(sample == Tube_code))%>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
beta_phylogenetic_nmds_post4 <- phylo_post4 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_post4, by = join_by(sample == Tube_code))%>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
ggarrange(beta_richness_nmds_post4, beta_neutral_nmds_post4, beta_phylogenetic_nmds_post4, ncol=2, nrow=2, common.legend = TRUE, legend="right")
Functional differences
WC from acclimation to FMT1
element_gift %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Hot_control") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 0 × 6
# ℹ 6 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>, Domain <chr>, Function <chr>, Element <chr>
What is the effect of FMT and temperature on the GM after 2 weeks caompared to acclimation?
CI vs CC
Alpha
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(type %in% c("Control","Treatment") & time_point %in% c("FMT2","Acclimation")) %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5, show.legend = FALSE) +
facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
# facet_nested(.~metric+type, labeller = labeller(metric = label_map, type = label_map))+
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
values=c("#4477AA", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control", "Treatment"),
values=c("#4477AA50","#76b18350")) +
theme_minimal() +
theme(axis.text.x=element_text(size=10))+
labs(x = "Time Point", y = "value")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Hot_control")
Richness
Modelq0GLMMNB <- MASS::glm.nb(richness ~ type*time_point, data = alpha_div_meta,trace=TRUE)
summary(Modelq0GLMMNB)
Call:
MASS::glm.nb(formula = richness ~ type * time_point, data = alpha_div_meta,
trace = TRUE, init.theta = 4.087439892, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.9512 0.1816 21.756 <2e-16 ***
typeTreatment -0.1372 0.2502 -0.548 0.584
time_pointFMT2 0.1206 0.2491 0.484 0.628
typeTreatment:time_pointFMT2 0.4149 0.3469 1.196 0.232
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(4.0874) family taken to be 1)
Null deviance: 43.336 on 34 degrees of freedom
Residual deviance: 37.889 on 31 degrees of freedom
AIC: 340.92
Number of Fisher Scoring iterations: 1
Theta: 4.09
Std. Err.: 1.07
2 x log-likelihood: -330.917
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ type)
$emmeans
type emmean SE df asymp.LCL asymp.UCL
Control 4.01 0.125 Inf 3.77 4.26
Treatment 4.08 0.121 Inf 3.85 4.32
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Control - Treatment -0.0702 0.173 Inf -0.405 0.6855
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Neutral
Modelq1n <- lm(formula = neutral ~ type*time_point, data = alpha_div_meta)
summary(Modelq1n)
Call:
lm(formula = neutral ~ type * time_point, data = alpha_div_meta)
Residuals:
Min 1Q Median 3Q Max
-20.377 -6.596 2.212 6.374 21.492
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 21.373 3.672 5.820 2.05e-06 ***
typeTreatment -3.974 5.047 -0.787 0.43704
time_pointFMT2 2.345 5.047 0.465 0.64552
typeTreatment:time_pointFMT2 20.138 7.032 2.864 0.00745 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.39 on 31 degrees of freedom
Multiple R-squared: 0.4388, Adjusted R-squared: 0.3845
F-statistic: 8.08 on 3 and 31 DF, p-value: 0.0004037
Phylogenetic
Model_phylo <- lm(formula = phylogenetic ~ type*time_point, data = alpha_div_meta)
summary(Model_phylo)
Call:
lm(formula = phylogenetic ~ type * time_point, data = alpha_div_meta)
Residuals:
Min 1Q Median 3Q Max
-3.6469 -0.6788 0.0109 0.4658 2.9735
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.22734 0.52818 9.897 4.1e-11 ***
typeTreatment 0.31260 0.72591 0.431 0.670
time_pointFMT2 -0.45260 0.72591 -0.623 0.538
typeTreatment:time_pointFMT2 0.02719 1.01138 0.027 0.979
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.494 on 31 degrees of freedom
Multiple R-squared: 0.03742, Adjusted R-squared: -0.05573
F-statistic: 0.4017 on 3 and 31 DF, p-value: 0.7527
Beta diversity
Number of samples used
samples_to_keep_post5 <- sample_metadata %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Hot_control") %>%
select(Tube_code) %>%
pull()
subset_meta_post5 <- sample_metadata %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Hot_control")
subset_meta_post5$type_time<-interaction(subset_meta_post5$type, subset_meta_post5$time_point)
length(samples_to_keep_post5)
[1] 35
Richness
richness_post5 <- as.matrix(beta_q0n$S)
richness_post5 <- as.dist(richness_post5[rownames(richness_post5) %in% samples_to_keep_post5,
colnames(richness_post5) %in% samples_to_keep_post5])
betadisper(richness_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.004294 0.0042942 0.4761 999 0.488
Residuals 33 0.297625 0.0090189
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Treatment
Control 0.486
Treatment 0.49501
adonis2(richness_post5 ~ type*time_point,
data = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))),
permutations = 999,
strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
=======
.table td.tinytable_css_qnbqx8co8bmcd7x3kkzt, .table th.tinytable_css_qnbqx8co8bmcd7x3kkzt { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
.table td.tinytable_css_fvxtb9e31nid36v4k606, .table th.tinytable_css_fvxtb9e31nid36v4k606 { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Acclimation vs FMT1 |
1 |
1.135334 |
3.400826 |
0.09606629 |
0.001 |
0.001 |
** |
Neutral
neutral_post3 <- as.matrix(beta_q1n$S)
neutral_post3 <- as.dist(neutral_post3[rownames(neutral_post3) %in% samples_to_keep_post3,
colnames(neutral_post3) %in% samples_to_keep_post3])
betadisper(neutral_post3, subset_meta_post3$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.019686 0.0196859 2.8068 999 0.109
Residuals 32 0.224435 0.0070136
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Treatment
Control 0.1
Treatment 0.10361
adonis2(neutral_post3 ~ type,
data = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))),
permutations = 999,
strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| type |
1 |
0.4361022 |
0.04050485 |
1.350872 |
1 |
| Residual |
32 |
10.3305628 |
0.95949515 |
NA |
NA |
| Total |
33 |
10.7666650 |
1.00000000 |
NA |
NA |
<<<<<<< HEAD
subset_meta_post3_arrange <- column_to_rownames(subset_meta_post3, "Tube_code")
subset_meta_post3_arrange<-subset_meta_post3_arrange[labels(neutral_post3),]
pairwise <- pairwise.adonis(neutral_post3, subset_meta_post3_arrange$type, perm=999)
pairwise%>%
tt()
subset_meta_post5_arrange <- column_to_rownames(subset_meta_post5, "Tube_code")
subset_meta_post5_arrange<-subset_meta_post5_arrange[labels(richness_post5),]
pairwise <- pairwise.adonis(richness_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
tt()
=======
.table td.tinytable_css_upmfaiopedcmahslh2sq, .table th.tinytable_css_upmfaiopedcmahslh2sq { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
.table td.tinytable_css_grdma35xmar66z9luemd, .table th.tinytable_css_grdma35xmar66z9luemd { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
<<<<<<< HEAD
| Control vs Treatment |
=======
Control.Acclimation vs Treatment.Acclimation |
1 |
0.3657243 |
1.123239 |
0.06966584 |
0.250 |
1.000 |
|
| Control.Acclimation vs Control.FMT2 |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
1 |
0.4361022 |
1.350872 |
0.04050485 |
0.143 |
0.143 |
|
pairwise <- pairwise.adonis(neutral_post3, subset_meta_post3_arrange$time_point, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
<<<<<<< HEAD
| Acclimation vs FMT1 |
=======
Control.Acclimation vs Treatment.FMT2 |
1 |
0.8072604 |
2.940901 |
0.16392161 |
0.003 |
0.018 |
. |
| Treatment.Acclimation vs Control.FMT2 |
1 |
0.5959230 |
1.984036 |
0.11032208 |
0.002 |
0.012 |
. |
| Treatment.Acclimation vs Treatment.FMT2 |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
1 |
1.683105 |
5.929324 |
0.1563256 |
0.001 |
0.001 |
** |
<<<<<<< HEAD
Phylogenetic
phylo_post3 <- as.matrix(beta_q1p$S)
phylo_post3 <- as.dist(phylo_post3[rownames(phylo_post3) %in% samples_to_keep_post3,
colnames(phylo_post3) %in% samples_to_keep_post3])
betadisper(phylo_post3, subset_meta_post3$type) %>% permutest(., pairwise = TRUE)
Neutral
neutral_post5 <- as.matrix(beta_q1n$S)
neutral_post5 <- as.dist(neutral_post5[rownames(neutral_post5) %in% samples_to_keep_post5,
colnames(neutral_post5) %in% samples_to_keep_post5])
betadisper(neutral_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
<<<<<<< HEAD
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04697 0.046974 2.8076 999 0.112
Residuals 32 0.53540 0.016731
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Treatment
Control 0.123
Treatment 0.10357
adonis2(phylo_post3 ~ type,
data = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))),
permutations = 999,
strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()
adonis2(neutral_post5 ~ type,
data = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))),
permutations = 999,
strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()
=======
.table td.tinytable_css_n02xo7zp855yrblo2als, .table th.tinytable_css_n02xo7zp855yrblo2als { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
.table td.tinytable_css_41s8p3wwr78e837y3ifn, .table th.tinytable_css_41s8p3wwr78e837y3ifn { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| type |
1 |
0.03710337 |
0.01720754 |
0.5602822 |
1 |
<<<<<<< HEAD
=======
0.642454 |
0.06322634 |
2.227293 |
1 |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| Residual |
32 |
2.11912470 |
0.98279246 |
NA |
NA |
| Total |
33 |
2.15622808 |
1.00000000 |
NA |
NA |
<<<<<<< HEAD
subset_meta_post3_arrange <- column_to_rownames(subset_meta_post3, "Tube_code")
subset_meta_post3_arrange<-subset_meta_post3_arrange[labels(phylo_post3),]
pairwise <- pairwise.adonis(phylo_post3, subset_meta_post3_arrange$type, perm=999)
pairwise%>%
tt()
subset_meta_post5_arrange <- column_to_rownames(subset_meta_post5, "Tube_code")
subset_meta_post5_arrange<-subset_meta_post5_arrange[labels(neutral_post5),]
pairwise <- pairwise.adonis(neutral_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
tt()
=======
.table td.tinytable_css_z3btk0dexo1og70pzmf6, .table th.tinytable_css_z3btk0dexo1og70pzmf6 { border-bottom: solid #d3d8dc 0.1em; }
.table td.tinytable_css_ao0irp265zjiku95dnnu, .table th.tinytable_css_ao0irp265zjiku95dnnu { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Control vs Treatment |
1 |
<<<<<<< HEAD
0.03710337 |
0.5602822 |
0.01720754 |
0.678 |
0.678 |
=======
0.2759858 |
0.9928976 |
0.06208366 |
0.450 |
1.000 |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
|
pairwise <- pairwise.adonis(phylo_post3, subset_meta_post3_arrange$time_point, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Acclimation vs FMT1 |
1 |
0.3325691 |
5.835636 |
0.1542365 |
0.001 |
0.001 |
** |
dbRDA
#Richness
cca_ord <- capscale(formula = richness_post3 ~ subset_meta_post3$time_point* subset_meta_post3$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post3, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post3$time_pointFMT1" = "FMT1",
"subset_meta_post3$typeTreatment" = "Cold-intervention",
"subset_meta_post3$time_pointFMT1:subset_meta_post3$typeTreatment" = "Interaction in Cold-intervention")
beta_richness_nmds_post3 <-CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Neutral
cca_ord <- capscale(formula = neutral_post3 ~ subset_meta_post3$time_point* subset_meta_post3$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post3, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post3$time_pointFMT1" = "FMT1",
"subset_meta_post3$typeTreatment" = "Cold-intervention",
"subset_meta_post3$time_pointFMT1:subset_meta_post3$typeTreatment" = "Interaction in Cold-intervention")
beta_neutral_nmds_post3 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Phylogenetic
cca_ord <- capscale(formula = phylo_post3 ~ subset_meta_post3$time_point* subset_meta_post3$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post3, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post3$time_pointFMT1" = "FMT1",
"subset_meta_post3$typeTreatment" = "Cold-intervention",
"subset_meta_post3$time_pointFMT1:subset_meta_post3$typeTreatment" = "Interaction in Cold-intervention")
beta_phylogenetic_nmds_post3 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
ggarrange(beta_richness_nmds_post3, beta_neutral_nmds_post3, beta_phylogenetic_nmds_post3, ncol=3, nrow=1, common.legend = TRUE, legend="right")
Functional differences
CC from acclimation to FMT1
element_gift %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Control") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 11 × 6
trait p_value p_adjust Domain Function Element
<chr> <dbl> <dbl> <chr> <chr> <chr>
1 B0204 0.00370 0.0478 Biosynthesis Amino acid biosynthesis_Serine Serine
2 B0205 0.000165 0.0117 Biosynthesis Amino acid biosynthesis_Threonine Threonine
3 B0703 0.00370 0.0478 Biosynthesis Vitamin biosynthesis_Niacin (B3) Niacin (B3)
4 B0805 0.00156 0.0370 Biosynthesis Aromatic compound biosynthesis_Indole-3-acetate Indole-3-acetate
5 D0102 0.000329 0.0156 Degradation Lipid degradation_Fatty acid Fatty acid
6 D0212 0.000987 0.0350 Degradation Polysaccharide degradation_Arabinan Arabinan
7 D0304 0.00247 0.0438 Degradation Sugar degradation_D-Arabinose D-Arabinose
8 D0306 0.00156 0.0370 Degradation Sugar degradation_D-Xylose D-Xylose
9 D0308 0.00247 0.0438 Degradation Sugar degradation_L-Rhamnose L-Rhamnose
10 D0601 0.00370 0.0478 Degradation Nitrogen compound degradation_Nitrate Nitrate
11 D0708 0.000165 0.0117 Degradation Alcohol degradation_Phytol Phytol
CI from acclimation to FMT1
element_gift %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Treatment") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 2 × 6
trait p_value p_adjust Domain Function Element
<chr> <dbl> <dbl> <chr> <chr> <chr>
1 B0710 0.0000823 0.0119 Biosynthesis Vitamin biosynthesis_Phylloquinone (K1) Phylloquinone (K1)
2 D0102 0.000576 0.0418 Degradation Lipid degradation_Fatty acid Fatty acid
CI vs WC
Alpha
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Control") %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
geom_boxplot(outlier.shape = NA) +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5) +
facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
values=c("#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
values=c("#d57d2c50","#76b18350")) +
theme_minimal() +
theme(axis.text.x=element_text(size=8))+
labs(x = "Time Point", y = "value")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Control")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(5.0429) ( log )
Formula: richness ~ type * time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
338.8 348.1 -163.4 326.8 29
Scaled residuals:
Min 1Q Median 3Q Max
-1.85185 -0.65278 0.06667 0.54823 1.97668
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 1.315e-11 3.627e-06
Number of obs: 35, groups: individual, 18
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.4697 0.1527 29.279 < 2e-16 ***
typeTreatment -0.6557 0.2186 -2.999 0.00271 **
time_pointFMT1 -0.4208 0.2174 -1.936 0.05292 .
typeTreatment:time_pointFMT1 0.5484 0.3146 1.743 0.08131 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) typTrt t_FMT1
typeTretmnt -0.698
tim_pntFMT1 -0.702 0.490
typTr:_FMT1 0.485 -0.695 -0.691
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ type)
$emmeans
type emmean SE df asymp.LCL asymp.UCL
Hot_control 4.26 0.109 Inf 4.05 4.47
Treatment 3.88 0.114 Inf 3.65 4.10
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Hot_control - Treatment 0.381 0.157 Inf 2.425 0.0153
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ time_point)
$emmeans
time_point emmean SE df asymp.LCL asymp.UCL
Acclimation 4.14 0.109 Inf 3.93 4.36
FMT1 4.00 0.113 Inf 3.77 4.22
Results are averaged over the levels of: type
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Acclimation - FMT1 0.147 0.157 Inf 0.932 0.3512
Results are averaged over the levels of: type
Results are given on the log (not the response) scale.
Neutral
Modelq1n <- lme(neutral ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta,)
summary(Modelq1n)
Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
250.0982 258.7021 -119.0491
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 5.122922 8.546652
Fixed effects: neutral ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 44.60136 3.321472 16 13.428192 0e+00
typeTreatment -27.20238 4.697271 16 -5.791103 0e+00
time_pointFMT1 -18.86154 4.028930 15 -4.681526 3e-04
typeTreatment:time_pointFMT1 26.15805 5.809237 15 4.502838 4e-04
Correlation:
(Intr) typTrt t_FMT1
typeTreatment -0.707
time_pointFMT1 -0.606 0.429
typeTreatment:time_pointFMT1 0.421 -0.595 -0.694
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.903393212 -0.678106330 0.002306491 0.652929042 1.401732098
Number of Observations: 35
Number of Groups: 18
emmeans::emmeans(Modelq1n, pairwise ~ type)
$emmeans
type emmean SE df lower.CL upper.CL
Hot_control 35.2 2.64 15 29.5 40.8
Treatment 21.0 2.70 15 15.3 26.8
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Hot_control - Treatment 14.1 3.78 15 3.739 0.0020
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
emmeans::emmeans(Modelq1n, pairwise ~ time_point)
$emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 31.0 2.35 16 26.0 36.0
FMT1 25.2 2.42 15 20.1 30.4
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - FMT1 5.78 2.9 15 1.991 0.0650
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Phylogenetic
Model_phylo <- lme(phylogenetic ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Model_phylo)
Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
124.0635 132.6674 -56.03174
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.0002455305 1.282332
Fixed effects: phylogenetic ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 6.522583 0.4274440 16 15.259503 0.0000
typeTreatment -0.982643 0.6044971 16 -1.625554 0.1236
time_pointFMT1 -2.247960 0.6044971 15 -3.718727 0.0021
typeTreatment:time_pointFMT1 0.898121 0.8681429 15 1.034531 0.3173
Correlation:
(Intr) typTrt t_FMT1
typeTreatment -0.707
time_pointFMT1 -0.707 0.500
typeTreatment:time_pointFMT1 0.492 -0.696 -0.696
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.40177643 -0.54641173 -0.06022755 0.42683546 2.21670618
Number of Observations: 35
Number of Groups: 18
emmeans::emmeans(Model_phylo, pairwise ~ type)
$emmeans
type emmean SE df lower.CL upper.CL
Hot_control 5.40 0.302 15 4.75 6.04
Treatment 4.87 0.312 15 4.20 5.53
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Hot_control - Treatment 0.534 0.434 15 1.229 0.2379
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
emmeans::emmeans(Model_phylo, pairwise ~ time_point)
$emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 6.03 0.302 16 5.39 6.67
FMT1 4.23 0.312 15 3.57 4.90
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - FMT1 1.8 0.434 15 4.144 0.0009
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Beta diversity
Number of samples used
samples_to_keep_post4 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Control") %>%
select(Tube_code) %>%
pull()
subset_meta_post4 <- sample_metadata %>%
filter(time_point=="FMT1"|time_point=="Acclimation") %>%
filter(type!="Control")
subset_meta_post4$time_point<-as.factor(subset_meta_post4$time_point)
subset_meta_post4$type<-as.factor(subset_meta_post4$type)
length(samples_to_keep_post4)
[1] 35
Richness
richness_post4 <- as.matrix(beta_q0n$S)
richness_post4 <- as.dist(richness_post4[rownames(richness_post4) %in% samples_to_keep_post4,
colnames(richness_post4) %in% samples_to_keep_post4])
betadisper(richness_post4, subset_meta_post4$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.076276 0.076276 19.518 999 0.001 ***
Residuals 33 0.128967 0.003908
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
Hot_control 0.001
Treatment 0.00010134
adonis2(richness_post4 ~ type*time_point,
data = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))),
permutations = 999,
strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| type |
1 |
1.1360754 |
0.09995918 |
4.050611 |
0.001 |
| time_point |
1 |
0.9251516 |
0.08140076 |
3.298575 |
0.001 |
| type:time_point |
1 |
0.6095926 |
0.05363586 |
2.173467 |
0.002 |
| Residual |
31 |
8.6945735 |
0.76500420 |
NA |
NA |
| Total |
34 |
11.3653932 |
1.00000000 |
NA |
NA |
subset_meta_post4_arrange <- column_to_rownames(subset_meta_post4, "Tube_code")
subset_meta_post4_arrange<-subset_meta_post4_arrange[labels(richness_post4),]
pairwise <- pairwise.adonis(richness_post4, subset_meta_post4_arrange$type, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Treatment vs Hot_control |
1 |
1.136075 |
3.665004 |
0.09995918 |
0.001 |
0.001 |
** |
pairwise <- pairwise.adonis(richness_post4, subset_meta_post4_arrange$time_point, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Acclimation vs FMT1 |
1 |
0.9366646 |
2.963922 |
0.08241375 |
0.001 |
0.001 |
** |
Neutral
neutral_post4 <- as.matrix(beta_q1n$S)
neutral_post4 <- as.dist(neutral_post4[rownames(neutral_post4) %in% samples_to_keep_post4,
colnames(neutral_post4) %in% samples_to_keep_post4])
betadisper(neutral_post4, subset_meta_post4$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.087335 0.087335 15.924 999 0.002 **
Residuals 33 0.180987 0.005484
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
Hot_control 0.001
Treatment 0.00034565
adonis2(neutral_post4 ~ type,
data = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))),
permutations = 999,
strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| type |
1 |
1.332543 |
0.1214599 |
4.562313 |
1 |
| Residual |
33 |
9.638511 |
0.8785401 |
NA |
NA |
| Total |
34 |
10.971054 |
1.0000000 |
NA |
NA |
subset_meta_post4_arrange <- column_to_rownames(subset_meta_post4, "Tube_code")
subset_meta_post4_arrange<-subset_meta_post4_arrange[labels(neutral_post4),]
pairwise <- pairwise.adonis(neutral_post4, subset_meta_post4_arrange$type, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Treatment vs Hot_control |
1 |
1.332543 |
4.562313 |
0.1214599 |
0.001 |
0.001 |
** |
pairwise <- pairwise.adonis(neutral_post4, subset_meta_post4_arrange$time_point, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Acclimation vs FMT1 |
1 |
1.291025 |
4.401211 |
0.1176756 |
0.001 |
0.001 |
** |
Phylogenetic
phylo_post4 <- as.matrix(beta_q1p$S)
phylo_post4 <- as.dist(phylo_post4[rownames(phylo_post4) %in% samples_to_keep_post4,
colnames(phylo_post4) %in% samples_to_keep_post4])
betadisper(phylo_post4, subset_meta_post4$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.07659 0.076588 5.0159 999 0.035 *
Residuals 33 0.50388 0.015269
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
Hot_control 0.036
Treatment 0.031971
adonis2(phylo_post4 ~ type,
data = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))),
permutations = 999,
strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| type |
1 |
0.2070779 |
0.09373244 |
3.413088 |
1 |
| Residual |
33 |
2.0021670 |
0.90626756 |
NA |
NA |
| Total |
34 |
2.2092449 |
1.00000000 |
NA |
NA |
subset_meta_post4_arrange <- column_to_rownames(subset_meta_post4, "Tube_code")
subset_meta_post4_arrange<-subset_meta_post4_arrange[labels(phylo_post4),]
pairwise <- pairwise.adonis(phylo_post4, subset_meta_post4_arrange$type, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Treatment vs Hot_control |
1 |
0.2070779 |
3.413088 |
0.09373244 |
0.012 |
0.012 |
. |
pairwise <- pairwise.adonis(phylo_post4, subset_meta_post4_arrange$time_point, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Acclimation vs FMT1 |
1 |
0.2895801 |
4.978027 |
0.1310765 |
0.002 |
0.002 |
* |
RDA
#Richness
cca_ord <- capscale(formula = richness_post4 ~ subset_meta_post4$time_point* subset_meta_post4$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post4, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post4$time_pointFMT1" = "FMT1",
"subset_meta_post4$typeTreatment" = "Cold-intervention",
"subset_meta_post4$time_pointFMT1:subset_meta_post4$typeTreatment" = "Interaction in Cold-intervention")
beta_richness_nmds_post4 <-CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Neutral
cca_ord <- capscale(formula = neutral_post4 ~ subset_meta_post4$time_point* subset_meta_post4$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post4, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post4$time_pointFMT1" = "FMT1",
"subset_meta_post4$typeTreatment" = "Cold-intervention",
"subset_meta_post4$time_pointFMT1:subset_meta_post4$typeTreatment" = "Interaction in Cold-intervention")
beta_neutral_nmds_post4 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Phylogenetic
cca_ord <- capscale(formula = phylo_post4 ~ subset_meta_post4$time_point* subset_meta_post4$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post4, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post4$time_pointFMT1" = "FMT1",
"subset_meta_post4$typeTreatment" = "Cold-intervention",
"subset_meta_post4$time_pointFMT1:subset_meta_post4$typeTreatment" = "Interaction in Cold-intervention")
beta_phylogenetic_nmds_post4 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
ggarrange(beta_richness_nmds_post4, beta_neutral_nmds_post4, beta_phylogenetic_nmds_post4, ncol=3, nrow=1, common.legend = TRUE, legend="right")
Functional differences
WC from acclimation to FMT1
element_gift %>%
filter(time_point %in% c("FMT1","Acclimation") & type == "Hot_control") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 0 × 6
# ℹ 6 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>, Domain <chr>, Function <chr>, Element <chr>
What is the effect of FMT and temperature on the GM after 2 weeks caompared to acclimation?
CI vs CC
Alpha
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(type %in% c("Control","Treatment") & time_point %in% c("FMT2","Acclimation")) %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5, show.legend = FALSE) +
facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
# facet_nested(.~metric+type, labeller = labeller(metric = label_map, type = label_map))+
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
values=c("#4477AA", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control", "Treatment"),
values=c("#4477AA50","#76b18350")) +
theme_minimal() +
theme(axis.text.x=element_text(size=10))+
labs(x = "Time Point", y = "value")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Hot_control")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(4.0875) ( log )
Formula: richness ~ type * time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
342.9 352.2 -165.5 330.9 29
Scaled residuals:
Min 1Q Median 3Q Max
-1.7595 -0.2361 0.2110 0.4983 1.1980
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 6.519e-12 2.553e-06
Number of obs: 35, groups: individual, 18
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.9512 0.1816 21.756 <2e-16 ***
typeTreatment -0.1372 0.2502 -0.548 0.583
time_pointFMT2 0.1206 0.2491 0.484 0.628
typeTreatment:time_pointFMT2 0.4149 0.3469 1.196 0.232
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) typTrt t_FMT2
typeTretmnt -0.726
tim_pntFMT2 -0.729 0.529
typTr:_FMT2 0.524 -0.721 -0.718
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ type)
$emmeans
type emmean SE df asymp.LCL asymp.UCL
Control 4.01 0.125 Inf 3.77 4.26
Treatment 4.08 0.121 Inf 3.85 4.32
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Control - Treatment -0.0702 0.173 Inf -0.405 0.6855
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ time_point)
$emmeans
time_point emmean SE df asymp.LCL asymp.UCL
Acclimation 3.88 0.125 Inf 3.64 4.13
FMT2 4.21 0.120 Inf 3.98 4.45
Results are averaged over the levels of: type
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Acclimation - FMT2 -0.328 0.173 Inf -1.892 0.0585
Results are averaged over the levels of: type
Results are given on the log (not the response) scale.
Neutral
Modelq1n <- lme(neutral ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Modelq1n)
Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
253.6069 262.2108 -120.8034
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 3.268167 9.858563
Fixed effects: neutral ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 21.455542 3.670059 16 5.846103 0.0000
typeTreatment -4.056558 5.045308 16 -0.804026 0.4332
time_pointFMT2 2.262002 4.804331 15 0.470826 0.6445
typeTreatment:time_pointFMT2 20.220508 6.684284 15 3.025082 0.0085
Correlation:
(Intr) typTrt t_FMT2
typeTreatment -0.727
time_pointFMT2 -0.697 0.507
typeTreatment:time_pointFMT2 0.501 -0.684 -0.719
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.7676771 -0.6504757 0.2922823 0.6462932 2.0120364
Number of Observations: 35
Number of Groups: 18
emmeans::emmeans(Modelq1n, pairwise ~ type)
$emmeans
type emmean SE df lower.CL upper.CL
Control 22.6 2.64 15 17.0 28.2
Treatment 28.6 2.57 15 23.2 34.1
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Control - Treatment -6.05 3.68 15 -1.645 0.1208
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
emmeans::emmeans(Modelq1n, pairwise ~ time_point)
$emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 19.4 2.52 16 14.1 24.8
FMT2 31.8 2.45 15 26.6 37.0
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - FMT2 -12.4 3.34 15 -3.702 0.0021
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Phylogenetic
Model_phylo <- lme(phylogenetic ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Model_phylo)
Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
133.508 142.1119 -60.75399
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.2974732 1.463985
Fixed effects: phylogenetic ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 5.233354 0.5281280 16 9.909252 0.0000
typeTreatment 0.306587 0.7258724 16 0.422370 0.6784
time_pointFMT2 -0.458617 0.7121982 15 -0.643946 0.5293
typeTreatment:time_pointFMT2 0.033201 0.9917181 15 0.033479 0.9737
Correlation:
(Intr) typTrt t_FMT2
typeTreatment -0.728
time_pointFMT2 -0.715 0.521
typeTreatment:time_pointFMT2 0.514 -0.705 -0.718
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.403145672 -0.437916442 -0.005178256 0.313089271 1.956678742
Number of Observations: 35
Number of Groups: 18
emmeans::emmeans(Model_phylo, pairwise ~ type)
$emmeans
type emmean SE df lower.CL upper.CL
Control 5.00 0.370 15 4.22 5.79
Treatment 5.33 0.359 15 4.56 6.09
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Control - Treatment -0.323 0.515 15 -0.627 0.5400
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
emmeans::emmeans(Model_phylo, pairwise ~ time_point)
$emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 5.39 0.363 16 4.62 6.16
FMT2 4.94 0.352 15 4.19 5.70
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - FMT2 0.442 0.496 15 0.891 0.3868
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Beta diversity
Number of samples used
samples_to_keep_post5 <- sample_metadata %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Hot_control") %>%
select(Tube_code) %>%
pull()
subset_meta_post5 <- sample_metadata %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Hot_control")
subset_meta_post5$time_point<-as.factor(subset_meta_post5$time_point)
subset_meta_post5$type<-as.factor(subset_meta_post5$type)
length(samples_to_keep_post5)
[1] 35
Richness
richness_post5 <- as.matrix(beta_q0n$S)
richness_post5 <- as.dist(richness_post5[rownames(richness_post5) %in% samples_to_keep_post5,
colnames(richness_post5) %in% samples_to_keep_post5])
betadisper(richness_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.004294 0.0042942 0.4761 999 0.52
Residuals 33 0.297625 0.0090189
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Treatment
Control 0.524
Treatment 0.49501
adonis2(richness_post5 ~ type*time_point,
data = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))),
permutations = 999,
strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| type |
1 |
0.5174241 |
0.04775345 |
1.797587 |
0.001 |
| time_point |
1 |
0.9068994 |
0.08369841 |
3.150666 |
0.001 |
| type:time_point |
1 |
0.4878440 |
0.04502349 |
1.694822 |
0.028 |
| Residual |
31 |
8.9231562 |
0.82352465 |
NA |
NA |
| Total |
34 |
10.8353238 |
1.00000000 |
NA |
NA |
subset_meta_post5_arrange <- column_to_rownames(subset_meta_post5, "Tube_code")
subset_meta_post5_arrange<-subset_meta_post5_arrange[labels(richness_post5),]
pairwise <- pairwise.adonis(richness_post5, subset_meta_post5_arrange$type, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Control vs Treatment |
1 |
0.5174241 |
1.65489 |
0.04775345 |
0.022 |
0.022 |
. |
pairwise <- pairwise.adonis(richness_post5, subset_meta_post5_arrange$time_point, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Acclimation vs FMT2 |
1 |
0.9000618 |
2.989558 |
0.08306737 |
0.001 |
0.001 |
** |
Neutral
neutral_post5 <- as.matrix(beta_q1n$S)
neutral_post5 <- as.dist(neutral_post5[rownames(neutral_post5) %in% samples_to_keep_post5,
colnames(neutral_post5) %in% samples_to_keep_post5])
betadisper(neutral_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.034255 0.034255 3.9748 999 0.063 .
Residuals 33 0.284399 0.008618
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Treatment
Control 0.058
Treatment 0.054506
adonis2(neutral_post5 ~ type*time_point,
data = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))),
permutations = 999,
strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| type |
1 |
0.6424540 |
0.06322634 |
2.589936 |
0.001 |
| time_point |
1 |
1.1894033 |
0.11705370 |
4.794863 |
0.001 |
| type:time_point |
1 |
0.6395258 |
0.06293816 |
2.578132 |
0.010 |
| Residual |
31 |
7.6897933 |
0.75678179 |
NA |
NA |
| Total |
34 |
10.1611764 |
1.00000000 |
NA |
NA |
subset_meta_post5_arrange <- column_to_rownames(subset_meta_post5, "Tube_code")
subset_meta_post5_arrange<-subset_meta_post5_arrange[labels(neutral_post5),]
pairwise <- pairwise.adonis(neutral_post5, subset_meta_post5_arrange$type, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Control vs Treatment |
1 |
0.642454 |
2.227293 |
0.06322634 |
0.007 |
0.007 |
* |
pairwise <- pairwise.adonis(neutral_post5, subset_meta_post5_arrange$time_point, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Acclimation vs FMT2 |
1 |
1.172649 |
4.305202 |
0.1154049 |
0.001 |
0.001 |
** |
Phylogenetic
phylo_post5 <- as.matrix(beta_q1p$S)
phylo_post5 <- as.dist(phylo_post5[rownames(phylo_post5) %in% samples_to_keep_post5,
colnames(phylo_post5) %in% samples_to_keep_post5])
betadisper(phylo_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.02670 0.026703 1.6775 999 0.185
Residuals 33 0.52532 0.015919
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Treatment
Control 0.205
Treatment 0.20425
adonis2(phylo_post5 ~ type*time_point,
data = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))),
permutations = 999,
strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| type |
1 |
0.03644585 |
0.01867131 |
0.7636208 |
0.003 |
| time_point |
1 |
0.35441680 |
0.18156873 |
7.4258113 |
0.002 |
| type:time_point |
1 |
0.08154954 |
0.04177806 |
1.7086422 |
0.204 |
| Residual |
31 |
1.47955831 |
0.75798190 |
NA |
NA |
| Total |
34 |
1.95197051 |
1.00000000 |
NA |
NA |
subset_meta_post5_arrange <- column_to_rownames(subset_meta_post5, "Tube_code")
subset_meta_post5_arrange<-subset_meta_post5_arrange[labels(phylo_post5),]
pairwise <- pairwise.adonis(phylo_post5, subset_meta_post5_arrange$type, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Control vs Treatment |
1 |
0.03644585 |
0.6278766 |
0.01867131 |
0.603 |
0.603 |
|
pairwise <- pairwise.adonis(phylo_post5, subset_meta_post5_arrange$time_point, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Acclimation vs FMT2 |
1 |
0.3562316 |
7.366895 |
0.1824984 |
0.001 |
0.001 |
** |
dbRDA
#Richness
cca_ord <- capscale(formula = richness_post5 ~ subset_meta_post5$time_point* subset_meta_post5$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post5, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post5$time_pointFMT2" = "FMT2",
"subset_meta_post5$typeTreatment" = "Cold-intervention",
"subset_meta_post5$time_pointFMT2:subset_meta_post5$typeTreatment" = "Interaction in Cold-intervention")
beta_richness_nmds_post5 <-CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Neutral
cca_ord <- capscale(formula = neutral_post5 ~ subset_meta_post5$time_point* subset_meta_post5$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post5, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post5$time_pointFMT2" = "FMT2",
"subset_meta_post5$typeTreatment" = "Cold-intervention",
"subset_meta_post5$time_pointFMT2:subset_meta_post5$typeTreatment" = "Interaction in Cold-intervention")
beta_neutral_nmds_post5 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Phylogenetic
cca_ord <- capscale(formula = phylo_post5 ~ subset_meta_post5$time_point* subset_meta_post5$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post5, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post5$time_pointFMT2" = "FMT2",
"subset_meta_post5$typeTreatment" = "Cold-intervention",
"subset_meta_post5$time_pointFMT2:subset_meta_post5$typeTreatment" = "Interaction in Cold-intervention")
beta_phylogenetic_nmds_post5 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
ggarrange(beta_richness_nmds_post5, beta_neutral_nmds_post5, beta_phylogenetic_nmds_post5, ncol=3, nrow=1, common.legend = TRUE, legend="right")
Functional differences
CC from acclimation to FMT2
element_gift %>%
filter(time_point %in% c("FMT2","Acclimation") & type == "Control") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 0 × 6
# ℹ 6 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>, Domain <chr>, Function <chr>, Element <chr>
CI from acclimation to FMT2
element_gift %>%
filter(time_point %in% c("FMT2","Acclimation") & type == "Treatment") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 17 × 6
trait p_value p_adjust Domain Function Element
<chr> <dbl> <dbl> <chr> <chr> <chr>
1 B0103 0.000787 0.0230 Biosynthesis Nucleic acid biosynthesis_UDP/UTP UDP/UTP
2 B0214 0.000494 0.0230 Biosynthesis Amino acid biosynthesis_Glutamate Glutamate
3 D0102 0.000782 0.0230 Degradation Lipid degradation_Fatty acid Fatty acid
4 D0205 0.00399 0.0448 Degradation Polysaccharide degradation_Pectin Pectin
5 D0210 0.00123 0.0234 Degradation Polysaccharide degradation_Beta-mannan Beta-mannan
6 D0308 0.00276 0.0366 Degradation Sugar degradation_L-Rhamnose L-Rhamnose
7 D0512 0.000494 0.0230 Degradation Amino acid degradation_Histidine Histidine
8 D0513 0.00564 0.0484 Degradation Amino acid degradation_Tryptophan Tryptophan
9 D0601 0.000494 0.0230 Degradation Nitrogen compound degradation_Nitrate Nitrate
10 D0613 0.00564 0.0484 Degradation Nitrogen compound degradation_Taurine Taurine
11 D0706 0.00564 0.0484 Degradation Alcohol degradation_Ethylene glycol Ethylene glycol
12 D0708 0.00185 0.0270 Degradation Alcohol degradation_Phytol Phytol
13 D0801 0.00145 0.0234 Degradation Xenobiotic degradation_Toluene Toluene
14 D0802 0.00145 0.0234 Degradation Xenobiotic degradation_Xylene Xylene
15 D0901 0.00564 0.0484 Degradation Antibiotic degradation_Penicillin Penicillin
16 D0902 0.00108 0.0234 Degradation Antibiotic degradation_Carbapenem Carbapenem
17 S0105 0.00399 0.0448 Structure Cellular structure_Lipopolysaccharide Lipopolysaccharide
CI vs WC
Alpha
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Control") %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
geom_boxplot(outlier.shape = NA) +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5) +
facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
values=c("#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
values=c("#d57d2c50","#76b18350")) +
theme_minimal() +
theme(axis.text.x=element_text(size=8))+
labs(x = "Time Point", y = "value")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Control")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ type*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
emmeans(Modelq0GLMMNB, pairwise ~ type)
emmeans(Modelq0GLMMNB, pairwise ~ time_point)
Neutral
Modelq1n <- lme(neutral ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Modelq1n)
Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
255.5098 264.3042 -121.7549
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 1.317597 9.38306
Fixed effects: neutral ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 44.60136 3.158373 16 14.121626 0.0000
typeTreatment -27.20238 4.466614 16 -6.090157 0.0000
time_pointFMT2 -12.40101 4.423217 16 -2.803618 0.0127
typeTreatment:time_pointFMT2 34.88352 6.255373 16 5.576569 0.0000
Correlation:
(Intr) typTrt t_FMT2
typeTreatment -0.707
time_pointFMT2 -0.700 0.495
typeTreatment:time_pointFMT2 0.495 -0.700 -0.707
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.1054551 -0.5575213 0.2345834 0.5912103 2.0213400
Number of Observations: 36
Number of Groups: 18
emmeans(Modelq1n, pairwise ~ type)
$emmeans
type emmean SE df lower.CL upper.CL
Hot_control 38.4 2.25 16 33.6 43.2
Treatment 28.6 2.25 16 23.9 33.4
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Hot_control - Treatment 9.76 3.19 16 3.061 0.0075
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
emmeans(Modelq1n, pairwise ~ time_point)
$emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 31 2.23 16 26.3 35.7
FMT2 36 2.23 16 31.3 40.8
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - FMT2 -5.04 3.13 16 -1.612 0.1266
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Phylogenetic
Model_phylo <- lme(phylogenetic ~ type*time_point,
random = ~ 1 | individual,
data = alpha_div_meta)
summary(Model_phylo)
Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
132.8692 141.6636 -60.43459
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.7604222 1.204042
Fixed effects: phylogenetic ~ type * time_point
Value Std.Error DF t-value p-value
(Intercept) 6.522583 0.4746882 16 13.740774 0.0000
typeTreatment -0.982643 0.6713105 16 -1.463768 0.1626
time_pointFMT2 -1.066614 0.5675911 16 -1.879195 0.0786
typeTreatment:time_pointFMT2 0.641199 0.8026950 16 0.798807 0.4361
Correlation:
(Intr) typTrt t_FMT2
typeTreatment -0.707
time_pointFMT2 -0.598 0.423
typeTreatment:time_pointFMT2 0.423 -0.598 -0.707
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.90560105 -0.54341938 -0.06251794 0.39498845 1.94321495
Number of Observations: 36
Number of Groups: 18
emmeans(Model_phylo, pairwise ~ type)
$emmeans
type emmean SE df lower.CL upper.CL
Hot_control 5.99 0.381 16 5.18 6.80
Treatment 5.33 0.381 16 4.52 6.13
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Hot_control - Treatment 0.662 0.538 16 1.230 0.2364
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
emmeans(Model_phylo, pairwise ~ time_point)
$emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 6.03 0.336 16 5.32 6.74
FMT2 5.29 0.336 16 4.57 6.00
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - FMT2 0.746 0.401 16 1.859 0.0816
Results are averaged over the levels of: type
Degrees-of-freedom method: containment
Beta diversity
Number of samples used
samples_to_keep_post6 <- sample_metadata %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Control") %>%
select(Tube_code) %>%
pull()
subset_meta_post6 <- sample_metadata %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Control")
subset_meta_post6$time_point<-as.factor(subset_meta_post6$time_point)
subset_meta_post6$type<-as.factor(subset_meta_post6$type)
length(samples_to_keep_post6)
[1] 36
Richness
richness_post6 <- as.matrix(beta_q0n$S)
richness_post6 <- as.dist(richness_post6[rownames(richness_post6) %in% samples_to_keep_post6,
colnames(richness_post6) %in% samples_to_keep_post6])
betadisper(richness_post6, subset_meta_post6$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.054351 0.054351 9.1422 999 0.004 **
Residuals 34 0.202132 0.005945
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
Hot_control 0.006
Treatment 0.0047271
adonis2(richness_post6 ~ type*time_point,
data = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))),
permutations = 999,
strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| type |
1 |
<<<<<<< HEAD
1.2570893 |
0.11659822 |
4.856166 |
0.001 |
=======
0.4748999 |
1.9609749 |
0.11561688 |
0.002 |
0.012 |
. |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| time_point |
1 |
0.6574224 |
0.06097760 |
2.539639 |
0.002 |
| type:time_point |
1 |
<<<<<<< HEAD
0.5831993 |
0.05409322 |
2.252913 |
0.004 |
| Residual |
32 |
8.2836650 |
0.76833097 |
NA |
NA |
| Total |
35 |
10.7813760 |
1.00000000 |
NA |
NA |
=======
0.6311089 |
2.4041625 |
0.13063146 |
0.002 |
0.012 |
. |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
subset_meta_post6_arrange <- column_to_rownames(subset_meta_post6, "Tube_code")
subset_meta_post6_arrange<-subset_meta_post6_arrange[labels(richness_post6),]
pairwise <- pairwise.adonis(richness_post6, subset_meta_post6_arrange$type, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Treatment vs Hot_control |
1 |
1.257089 |
4.487584 |
0.1165982 |
0.001 |
0.001 |
** |
pairwise <- pairwise.adonis(richness_post6, subset_meta_post6_arrange$time_point, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Acclimation vs FMT2 |
1 |
0.6574224 |
2.207869 |
0.0609776 |
0.003 |
0.003 |
* |
<<<<<<< HEAD
Neutral
neutral_post6 <- as.matrix(beta_q1n$S)
neutral_post6 <- as.dist(neutral_post6[rownames(neutral_post6) %in% samples_to_keep_post6,
colnames(neutral_post6) %in% samples_to_keep_post6])
betadisper(neutral_post6, subset_meta_post6$type) %>% permutest(., pairwise = TRUE)
=======
Phylogenetic
phylo_post5 <- as.matrix(beta_q1p$S)
phylo_post5 <- as.dist(phylo_post5[rownames(phylo_post5) %in% samples_to_keep_post5,
colnames(phylo_post5) %in% samples_to_keep_post5])
betadisper(phylo_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
<<<<<<< HEAD
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.078328 0.078328 12.006 999 0.001 ***
Residuals 34 0.221814 0.006524
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
Hot_control 0.005
Treatment 0.0014541
adonis2(neutral_post6 ~ type,
data = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))),
permutations = 999,
strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()
=======
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.02670 0.026703 1.6775 999 0.211
Residuals 33 0.52532 0.015919
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Treatment
Control 0.233
Treatment 0.20425
adonis2(phylo_post5 ~ type,
data = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))),
permutations = 999,
strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
=======
.table td.tinytable_css_9qk78kjt0i9lfnp4ejpp, .table th.tinytable_css_9qk78kjt0i9lfnp4ejpp { border-bottom: solid #d3d8dc 0.1em; }
.table td.tinytable_css_25z8w88t8nonj4rkmnqn, .table th.tinytable_css_25z8w88t8nonj4rkmnqn { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| type |
1 |
1.325563 |
0.1261942 |
4.910249 |
1 |
<<<<<<< HEAD
=======
0.03644585 |
0.01867131 |
0.6278766 |
1 |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| Residual |
34 |
9.178585 |
0.8738058 |
NA |
NA |
| Total |
35 |
10.504148 |
1.0000000 |
NA |
NA |
<<<<<<< HEAD
subset_meta_post6_arrange <- column_to_rownames(subset_meta_post6, "Tube_code")
subset_meta_post6_arrange<-subset_meta_post6_arrange[labels(neutral_post6),]
pairwise <- pairwise.adonis(neutral_post6, subset_meta_post6_arrange$type, perm=999)
pairwise%>%
tt()
NMDS
beta_richness_nmds_post5 <- richness_post5 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_post5, by = join_by(sample == Tube_code)) %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
beta_neutral_nmds_post5 <- neutral_post5 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_post5, by = join_by(sample == Tube_code))%>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
beta_phylogenetic_nmds_post5 <- phylo_post5 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_post5, by = join_by(sample == Tube_code))%>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Treatment"),
labels=c("Cold-Cold", "Cold-Hot"),
values=c("#4477AA","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
ggarrange(beta_richness_nmds_post5, beta_neutral_nmds_post5, beta_phylogenetic_nmds_post5, ncol=2, nrow=2, common.legend = TRUE, legend="right")
Functional differences
CC from acclimation to FMT2
element_gift %>%
filter(time_point %in% c("FMT2","Acclimation") & type == "Control") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 0 × 6
# ℹ 6 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>, Domain <chr>, Function <chr>, Element <chr>
CI from acclimation to FMT2
element_gift %>%
filter(time_point %in% c("FMT2","Acclimation") & type == "Treatment") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
# A tibble: 17 × 6
trait p_value p_adjust Domain Function Element
<chr> <dbl> <dbl> <chr> <chr> <chr>
1 B0103 0.000787 0.0230 Biosynthesis Nucleic acid biosynthesis_UDP/UTP UDP/UTP
2 B0214 0.000494 0.0230 Biosynthesis Amino acid biosynthesis_Glutamate Glutamate
3 D0102 0.000782 0.0230 Degradation Lipid degradation_Fatty acid Fatty acid
4 D0205 0.00399 0.0448 Degradation Polysaccharide degradation_Pectin Pectin
5 D0210 0.00123 0.0234 Degradation Polysaccharide degradation_Beta-mannan Beta-mannan
6 D0308 0.00276 0.0366 Degradation Sugar degradation_L-Rhamnose L-Rhamnose
7 D0512 0.000494 0.0230 Degradation Amino acid degradation_Histidine Histidine
8 D0513 0.00564 0.0484 Degradation Amino acid degradation_Tryptophan Tryptophan
9 D0601 0.000494 0.0230 Degradation Nitrogen compound degradation_Nitrate Nitrate
10 D0613 0.00564 0.0484 Degradation Nitrogen compound degradation_Taurine Taurine
11 D0706 0.00564 0.0484 Degradation Alcohol degradation_Ethylene glycol Ethylene glycol
12 D0708 0.00185 0.0270 Degradation Alcohol degradation_Phytol Phytol
13 D0801 0.00145 0.0234 Degradation Xenobiotic degradation_Toluene Toluene
14 D0802 0.00145 0.0234 Degradation Xenobiotic degradation_Xylene Xylene
15 D0901 0.00564 0.0484 Degradation Antibiotic degradation_Penicillin Penicillin
16 D0902 0.00108 0.0234 Degradation Antibiotic degradation_Carbapenem Carbapenem
17 S0105 0.00399 0.0448 Structure Cellular structure_Lipopolysaccharide Lipopolysaccharide
CI vs WC
Alpha
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Control") %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic"))) %>%
ggplot(aes(y = value, x = time_point, color = type, fill = type)) +
geom_boxplot(outlier.shape = NA) +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.5) +
facet_grid(metric ~ type, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
values=c("#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
values=c("#d57d2c50","#76b18350")) +
theme_minimal() +
theme(axis.text.x=element_text(size=8))+
labs(x = "Time Point", y = "value")
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Control")
Richness
Modelq0GLMMNB <- glm.nb(richness ~ type*time_point, data = alpha_div_meta,trace=TRUE)
summary(Modelq0GLMMNB)
emmeans(Modelq0GLMMNB, pairwise ~ type)
Neutral
Modelq1n <- lm(formula = neutral ~ type*time_point, data = alpha_div_meta)
summary(Modelq1n)
Call:
lm(formula = neutral ~ type * time_point, data = alpha_div_meta)
Residuals:
Min 1Q Median 3Q Max
-20.377 -5.358 2.158 5.663 19.246
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 44.601 3.158 14.122 2.68e-15 ***
typeTreatment -27.202 4.467 -6.090 8.36e-07 ***
time_pointFMT2 -12.401 4.467 -2.776 0.00911 **
typeTreatment:time_pointFMT2 34.884 6.317 5.522 4.34e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.475 on 32 degrees of freedom
Multiple R-squared: 0.571, Adjusted R-squared: 0.5308
F-statistic: 14.2 on 3 and 32 DF, p-value: 4.69e-06
Phylogenetic
Model_phylo <- lm(formula = phylogenetic ~ type*time_point, data = alpha_div_meta)
summary(Model_phylo)
Call:
lm(formula = phylogenetic ~ type * time_point, data = alpha_div_meta)
Residuals:
Min 1Q Median 3Q Max
-3.0799 -0.7676 -0.0971 0.5549 2.9735
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.5226 0.4747 13.741 5.7e-15 ***
typeTreatment -0.9826 0.6713 -1.464 0.153
time_pointFMT2 -1.0666 0.6713 -1.589 0.122
typeTreatment:time_pointFMT2 0.6412 0.9494 0.675 0.504
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.424 on 32 degrees of freedom
Multiple R-squared: 0.1321, Adjusted R-squared: 0.05075
F-statistic: 1.624 on 3 and 32 DF, p-value: 0.2033
Beta diversity
Number of samples used
samples_to_keep_post6 <- sample_metadata %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Control") %>%
select(Tube_code) %>%
pull()
subset_meta_post6 <- sample_metadata %>%
filter(time_point=="FMT2"|time_point=="Acclimation") %>%
filter(type!="Control")
subset_meta_post6$type_time<-interaction(subset_meta_post6$type, subset_meta_post6$time_point)
length(samples_to_keep_post6)
[1] 36
Richness
richness_post6 <- as.matrix(beta_q0n$S)
richness_post6 <- as.dist(richness_post6[rownames(richness_post6) %in% samples_to_keep_post6,
colnames(richness_post6) %in% samples_to_keep_post6])
betadisper(richness_post6, subset_meta_post6$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.054351 0.054351 9.1422 999 0.005 **
Residuals 34 0.202132 0.005945
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
Hot_control 0.005
Treatment 0.0047271
adonis2(richness_post6 ~ type*time_point,
data = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))),
permutations = 999,
strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
=======
.table td.tinytable_css_f0zf6n7lz8158paups9l, .table th.tinytable_css_f0zf6n7lz8158paups9l { border-bottom: solid #d3d8dc 0.1em; }
.table td.tinytable_css_008xqo3uvzkbpi6zctrt, .table th.tinytable_css_008xqo3uvzkbpi6zctrt { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Treatment vs Hot_control |
1 |
1.325563 |
4.910249 |
0.1261942 |
0.001 |
0.001 |
** |
<<<<<<< HEAD
pairwise <- pairwise.adonis(neutral_post6, subset_meta_post6_arrange$time_point, perm=999)
pairwise%>%
tt()
subset_meta_post6_arrange <- column_to_rownames(subset_meta_post6, "Tube_code")
subset_meta_post6_arrange<-subset_meta_post6_arrange[labels(richness_post6),]
pairwise <- pairwise.adonis(richness_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
tt()
=======
.table td.tinytable_css_ruz9z9h43gtq848ow32t, .table th.tinytable_css_ruz9z9h43gtq848ow32t { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
.table td.tinytable_css_ob11bhwyfi3ccm51mjew, .table th.tinytable_css_ob11bhwyfi3ccm51mjew { border-bottom: solid #d3d8dc 0.1em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
<<<<<<< HEAD
| Acclimation vs FMT2 |
=======
Treatment.Acclimation vs Hot_control.Acclimation |
1 |
1.3606630 |
5.087152 |
0.24124415 |
0.001 |
0.006 |
* |
| Treatment.Acclimation vs Treatment.FMT2 |
1 |
0.9130048 |
3.195028 |
0.16645080 |
0.001 |
0.006 |
* |
| Treatment.Acclimation vs Hot_control.FMT2 |
1 |
1.2747787 |
4.275366 |
0.21086503 |
0.001 |
0.006 |
* |
| Hot_control.Acclimation vs Treatment.FMT2 |
1 |
0.6397330 |
2.913695 |
0.15405213 |
0.001 |
0.006 |
* |
| Hot_control.Acclimation vs Hot_control.FMT2 |
1 |
0.3276169 |
1.412318 |
0.08111028 |
0.042 |
0.252 |
|
| Treatment.FMT2 vs Hot_control.FMT2 |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
1 |
0.9072453 |
3.214197 |
0.0863702 |
0.002 |
0.002 |
* |
<<<<<<< HEAD
Phylogenetic
phylo_post6 <- as.matrix(beta_q1p$S)
phylo_post6 <- as.dist(phylo_post6[rownames(phylo_post6) %in% samples_to_keep_post6,
colnames(phylo_post6) %in% samples_to_keep_post6])
betadisper(phylo_post6, subset_meta_post6$type) %>% permutest(., pairwise = TRUE)
Neutral
neutral_post6 <- as.matrix(beta_q1n$S)
neutral_post6 <- as.dist(neutral_post6[rownames(neutral_post6) %in% samples_to_keep_post6,
colnames(neutral_post6) %in% samples_to_keep_post6])
betadisper(neutral_post6, subset_meta_post6$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
<<<<<<< HEAD
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04938 0.049381 3.8729 999 0.056 .
Residuals 34 0.43351 0.012750
=======
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.078328 0.078328 12.006 999 0.002 **
Residuals 34 0.221814 0.006524
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
<<<<<<< HEAD
Hot_control 0.04
Treatment 0.057271
adonis2(phylo_post6 ~ type,
data = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))),
permutations = 999,
strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()
adonis2(neutral_post6 ~ type,
data = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))),
permutations = 999,
strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()
=======
.table td.tinytable_css_ov8p6edmajbwrv2eq9zh, .table th.tinytable_css_ov8p6edmajbwrv2eq9zh { border-bottom: solid #d3d8dc 0.1em; }
.table td.tinytable_css_22rqyvvcvazqw37ciq4q, .table th.tinytable_css_22rqyvvcvazqw37ciq4q { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| type |
1 |
0.1600877 |
0.08395907 |
3.116245 |
1 |
<<<<<<< HEAD
=======
1.325563 |
0.1261942 |
4.910249 |
1 |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| Residual |
34 |
1.7466474 |
0.91604093 |
NA |
NA |
| Total |
35 |
1.9067351 |
1.00000000 |
NA |
NA |
<<<<<<< HEAD
subset_meta_post6_arrange <- column_to_rownames(subset_meta_post6, "Tube_code")
subset_meta_post6_arrange<-subset_meta_post6_arrange[labels(phylo_post6),]
pairwise <- pairwise.adonis(phylo_post6, subset_meta_post6_arrange$type, perm=999)
pairwise%>%
tt()
subset_meta_post6_arrange <- column_to_rownames(subset_meta_post6, "Tube_code")
subset_meta_post6_arrange<-subset_meta_post6_arrange[labels(neutral_post6),]
pairwise <- pairwise.adonis(neutral_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
tt()
=======
.table td.tinytable_css_vj532nvb8indd1jitiw0, .table th.tinytable_css_vj532nvb8indd1jitiw0 { border-bottom: solid #d3d8dc 0.1em; }
.table td.tinytable_css_9qbyiug1baoth4sa8w1r, .table th.tinytable_css_9qbyiug1baoth4sa8w1r { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Treatment vs Hot_control |
1 |
<<<<<<< HEAD
0.1600877 |
3.116245 |
0.08395907 |
0.016 |
0.016 |
. |
=======
1.6125755 |
5.982598 |
0.27215155 |
0.001 |
0.006 |
* |
| Hot_control.Acclimation vs Treatment.FMT2 |
1 |
0.6202327 |
3.151987 |
0.16457754 |
0.001 |
0.006 |
* |
| Hot_control.Acclimation vs Hot_control.FMT2 |
1 |
0.3634438 |
1.708339 |
0.09647087 |
0.026 |
0.156 |
|
| Treatment.FMT2 vs Hot_control.FMT2 |
1 |
0.5010202 |
2.206532 |
0.12119453 |
0.001 |
0.006 |
* |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
<<<<<<< HEAD
pairwise <- pairwise.adonis(phylo_post6, subset_meta_post6_arrange$time_point, perm=999)
pairwise%>%
tt()
Phylogenetic
phylo_post6 <- as.matrix(beta_q1p$S)
phylo_post6 <- as.dist(phylo_post6[rownames(phylo_post6) %in% samples_to_keep_post6,
colnames(phylo_post6) %in% samples_to_keep_post6])
betadisper(phylo_post6, subset_meta_post6$type) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04938 0.049381 3.8729 999 0.048 *
Residuals 34 0.43351 0.012750
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Hot_control Treatment
Hot_control 0.034
Treatment 0.057271
adonis2(phylo_post6 ~ type,
data = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))),
permutations = 999,
strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()
=======
.table td.tinytable_css_nfn92de6sdn79asszwrh, .table th.tinytable_css_nfn92de6sdn79asszwrh { border-bottom: solid #d3d8dc 0.1em; }
.table td.tinytable_css_7xwp1sx2dsfixwp5fgpa, .table th.tinytable_css_7xwp1sx2dsfixwp5fgpa { border-top: solid #d3d8dc 0.1em; border-bottom: solid #d3d8dc 0.05em; }
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Acclimation vs FMT2 |
1 |
<<<<<<< HEAD
0.2882491 |
6.055333 |
0.1511742 |
0.001 |
0.001 |
** |
=======
0.1600877 |
0.08395907 |
3.116245 |
1 |
| Residual |
34 |
1.7466474 |
0.91604093 |
NA |
NA |
| Total |
35 |
1.9067351 |
1.00000000 |
NA |
NA |
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
<<<<<<< HEAD
dbRDA
#Richness
cca_ord <- capscale(formula = richness_post6 ~ subset_meta_post6$time_point* subset_meta_post6$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post6, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post6$time_pointFMT2" = "FMT2",
"subset_meta_post6$typeTreatment" = "Cold-intervention",
"subset_meta_post6$time_pointFMT2:subset_meta_post6$typeTreatment" = "Interaction in Cold-intervention")
beta_richness_nmds_post6 <-CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Neutral
cca_ord <- capscale(formula = neutral_post6 ~ subset_meta_post6$time_point* subset_meta_post6$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post6, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post6$time_pointFMT2" = "FMT2",
"subset_meta_post6$typeTreatment" = "Cold-intervention",
"subset_meta_post6$time_pointFMT2:subset_meta_post6$typeTreatment" = "Interaction in Cold-intervention")
beta_neutral_nmds_post6 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Phylogenetic
cca_ord <- capscale(formula = phylo_post6 ~ subset_meta_post6$time_point* subset_meta_post6$type)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta_post6, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta_post6$time_pointFMT2" = "FMT2",
"subset_meta_post6$typeTreatment" = "Cold-intervention",
"subset_meta_post6$time_pointFMT2:subset_meta_post6$typeTreatment" = "Interaction in Cold-intervention")
beta_phylogenetic_nmds_post6 <- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=type,shape = time_point)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
ggarrange(beta_richness_nmds_post6, beta_neutral_nmds_post6, beta_phylogenetic_nmds_post6, ncol=3, nrow=1, common.legend = TRUE, legend="right")
=======
NMDS
beta_richness_nmds_post6 <- richness_post6 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_post6, by = join_by(sample == Tube_code)) %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
beta_neutral_nmds_post6 <- neutral_post6 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_post6, by = join_by(sample == Tube_code))%>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
beta_phylogenetic_nmds_post6 <- phylo_post6 %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace = FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_post6, by = join_by(sample == Tube_code))%>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Hot_control", "Treatment"),
labels=c("Warm-control", "Cold-intervention"),
values=c("#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
ggarrange(beta_richness_nmds_post6, beta_neutral_nmds_post6, beta_phylogenetic_nmds_post6, ncol=2, nrow=2, common.legend = TRUE, legend="right")
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
Functional differences
WC from acclimation to FMT2
<<<<<<< HEAD
element_gift %>%
filter(time_point %in% c("FMT2","Acclimation") & type == "Hot_control") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
=======
element_gift %>%
filter(time_point %in% c("FMT2","Acclimation") & type == "Hot_control") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 10)], by = "Tube_code") %>%
pivot_longer(-c(Tube_code, time_point),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
remove_rownames() %>%
left_join(., uniqueGIFT_db, by=join_by("trait"=="Code_element"))
>>>>>>> d08d3ecaca65f0411fd9dfdf7dedf41ce523b6b2
# A tibble: 0 × 6
# ℹ 6 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>, Domain <chr>, Function <chr>, Element <chr>